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Abstract 


Hamiltonian structure is pursued and uncovered in collisional and collisionless gyrokinetic 
theory. A new Hamiltonian formulation of collisionless electromagnetic theory is presented that 
is ideally suited to implementation on modem supercomputers. The method used to uncover this 
structure is described in detail and applied to a number of examples, where several well-known 
plasma models are endowed with a Hamiltonian structure for the first time. The first energy- and 
momentum-conserving formulation of full-F collisional gyrokinetics is presented. In an effort to 
understand the theoretical underpinnings of this result at a deeper level, a stochastic Hamiltonian 
modeling approach is presented and applied to pitch angle scattering. Interestingly, the collision 
operator produced by the Hamiltonian approach is equal to the Lorentz operator plus higher-order 
terms, but does not exactly conserve energy. Conversely, the classical Lorentz collision operator is 
provably not Hamiltonian in the stochastic sense. 
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Chapter 1 
Introduction 


On microscopic scales, the physics of plasmas is Hamiltonian in nature. Neglecting quantum, rela¬ 
tivistic, and radiative effects for simplicity’s sake, the microscopic description of a plasma consists 
of a separate instance of the Lorentz force law for each plasma particle along with Maxwell’s 
equations to couple everything together. This system of equations can be derived from a varia¬ 
tional principle, which in turn can be used to derive a microscopic Hamiltonian functional and 
Poisson bracket. These equations governing microscopic plasma physics have been called the 
Klimontovich-Maxwell system, and their Lagrangian formulation is described for instance in 
Qin et al. (2014). 

With Hamiltonian structure ingrained so deeply in the foundation of the subject, plasma theory 
ought to be some grand exercise in the broader theory of Hamiltonian systems. And when viewed 
from a great distance, it is! However, for those in the trenches, studying plasma theory on a day- 
to-day basis, things seem different. The equations governing microscopic plasma theory are so 
hopelessly complicated that reduced plasma models are typically preferable to the Klimontovich- 
Maxwell model. These reduced models are obtained by carefully and cleverly applying Occam’s 
razor in order to tame the mathematical morass presented by the microscopic equations of motion. 
Sometimes, for instance in the case of ideal magnetohydrodynamics, the reduced model is prov- 


1 


Chapter 1. Introduction 


2 


ably Hamiltonian in nature Morrison and Greene (1980) However, often times the eonneetion 
between the reduced model and Hamiltonian mechanics is hazy at best. The most striking examples 
of this divorce from Hamiltonian mechanics arise when collisions must be accounted for within 
the confines of a continuum model. For instance the Landau collision operator probably cannot 
be derived from a conventional variational principle. The same is true of the more-fundamental 
Balescu-Lenard collision operator. More generally, there is a common feeling amongst plasma 
physicists that “dissipative dynamics are not Hamiltonian.” Even when collisions are neglected, 
there are reduced models that either fail to be Hamiltonian in nature, or so far have resisted at¬ 
tempts to be cast in Hamiltonian form. 

This thesis is devoted to revealing some new connections between Hamiltonian mechanics 
and the particularly interesting reduced plasma model known as gyrokinetics. I believe that it 
illustrates two general points in connection with the “Hamiltonian dichotomy” between reduced 
plasma models and the microscopic plasma model just described. First, Hamiltonian mechanics 
can be surprisingly useful in the study of reduced plasma models; the benefits of exploiting the 
Hamiltonian formalism to formulate and study reduced models can be unpredictable. Second, the 
connection between collisional plasma models and Hamiltonian mechanics is surprisingly deep; 
while collision operators do not fit within the traditional Hamiltonian framework, they may very 
well fit within a stochastic Hamiltonian framework Lazaro-Camf and Ortega (2008). The first 
point is covered roughly by Chapters 2 through 4, while the second point is discussed in Chapter 
5. 

In Chapter 2,1 present the results of an attempt to cast collisionless electromagnetic gyrokinet¬ 
ics in Hamiltonian form, a theory that already enjoys several Lagrangian formulations. The earliest 
of the Lagrangian formulations are given in Sugama (2000); Brizard (2000b,a), while more recent 
additions can be found in Pfirsch and Correa-Restrepo (2004); Squire et al. (2013). Given the typ- 

*In the case of ideal MHD, as well as many other Hamiltonian reduced models, it is still unclear how the reduced 
model’s Hamiltonian structure is related to the microscopic Hamiltonian structure. 
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ical intimate relationship between the Hamiltonian and Lagrangian formalisms, it is tempting to 
believe that the Hamiltonian structure of electromagnetic gyrokinetics should be straightforward 
to obtain. This was not the case. The usual Legendre transform technique fails when applied to 
gyrokinetic Lagrangians (one problem is these Lagrangians are linear in the phase space velocity 
of plasma species). Nevertheless, the existence of Lagrangian formulations of electromagnetic 
gyrokinetics suggested that a Hamiltonian formulation should exist; the riddle was how to find 
it. In the process of searching for this Hamiltonian structure, and with guidance from Professor 
Morrison, I found it technically convenient to slightly reformulate electromagnetic gyrokinetics in 
order to work with a manifestly gauge-invarianC theory along the lines of Morrison (2013). Sur¬ 
prisingly, this reformulation turned out to have several features that make it especially well-suited 
to simulation on modem supercomputers. This is the first example in the thesis of a surprising 
consequence of pursuing Hamiltonian structure in reduced plasma models. 

In chapter 31 give an extended account of the theoretical machine used to derive the gyrokinetic 
Poisson bracket from Chapter 2. This machine, which eats (possibly degenerate) Lagrangians and 
spits out Poisson brackets, is very closely related to the Peierls bracket formalism Peierls (1952), 
as well as the Dirac constraint formalism used in Squire et al. (2013). Nevertheless, several of the 
examples worked out in this chapter (besides electromagnetic gyrokinetics) are new. The first ex¬ 
ample that contains a new result concerns the Vlasov-Darwin system. This system had previously 
been cast in Hamiltonian form in Krause et al. (2007) using position-canonical momentum coor¬ 
dinates on the single-particle phase space. The novelty of the example in this Chapter is that the 
derivation of the bracket is done using position-velocity coordinates on the single-particle phase 
space (which leads to a different expression for the bracket.) The second novel example is con¬ 
cerned with deriving a bracket for the (quasi) neutral Vlasov system introduced by C. Tronci and 
E. Camporeale in Tronci and Camporeale (2015). Tronci and Camporeale provide a Lagrangian 
formulation of this reduced model, but stop short of passing to the Hamiltonian side. Thus, this 

^Previous work on Lagrangian electromagnetic gyrokinetics was usually done in the Coulomb gauge. 
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example presents a Poisson bracket and Hamiltonian functional for the neutral Vlasov model for 
the first time. The final novel example derives a Poisson bracket for force-free electrodynamics 
Gralla and Jacobsen (2014). Previously, a canonical bracket for this system was given in terms of 
an Euler potential representation of the electromagnetic field. The non-canonical bracket presented 
here does not use the Euler potential representation, and so is more general and potentially more 
useful for future stability analyses. 

Chapter 4 presents a result obtained in conjunction with Professor Brizard on an energy and 
momentum conserving nonlinear collision operator for full-f gyrokinetics. This work serves as 
the second example of a surprising consequence of applying the Hamiltonian formalism. The 
key insight that lead to this collision operator was a peculiar way of expressing the particle-space 
Eandau operator in terms of single-particle Poisson brackets. A deep reason as to why the Poisson 
bracket representation is as useful as it appears to be is still missing, and this is why the result is 
somewhat surprising. 

Einally, Chapter 5 contains the beginnings of a search for Hamiltonian structure underlying 
reduced collisional models at a deeper level than discussed in Chapter 4. It begins by describing 
generally how stochastic acceleration problems can be formulated in terms of stochastic Hamil¬ 
tonian mechanics Lazaro-Cami and Ortega (2008). Where ordinary Hamiltonian mechanics is 
concerned with one-parameter subgroups of the phase space symplectomorphism group, stochas¬ 
tic Hamiltonian mechanics is concerned with Brownian motion on the symplectomorphism group. 
The remainder of the chapter is then devoted to applying these ideas to the pitch angle scattering 
problem, which can be formulated as an example of stochastic acceleration. A surprise here is that 
there is a tension between energy conservation and the stochastic Hamiltonian formalism. The col¬ 
lision operator produced using the Hamiltonian approach is equal to the Eorentz collision operator 
plus higher-order terms, but does not exactly conserve kinetic energy. Conversely, the classical 
Eorentz operator is provably outside the realm of stochastic Hamiltonian mechanics. As I discuss 
at the end of the chapter, it seems likely that a way to overcome this problem is to slighly relax the 
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Markov approximation. However, this realization came too late in my stay at Princeton, and so I 
was not able to explore it more fully in this thesis. 




Chapter 2 

Hamiltonian formulation of the gyrokinetic 
Vlasov-Maxwell equations 

2.1 Introduction 

Electromagnetic gyrokinetic theory (EMGT) is a model used to describe the turbulent transport of 
particles and heat induced by fluctuating electric and magnetic fields in strongly magnetized plas¬ 
mas. EMGT is, in many ways, a more utilitarian tool than the more-fundamental Vlasov-Maxwell 
kinetic theory (VMKT). However, VMKT enjoys two important advantages over existing formu¬ 
lations of EMGT. (I) When simulated on a computer, the VMKT field solve is local; advancing 
the electromagnetic field in time at a given grid point only requires communication with nearby 
grid points Bowers et al. (2009). (II) There is an energy principle for assessing the stability of 
Vlasov-Maxwell equilibria Morrison and Pfirsch (1989) (also see Kruskal and Oberman (1958); 
Holm et al. (1985); Morrison (1998); Andreussi et al. (2012, 2013) for similar energy principles in 
other contexts). In contrast, modem EMGT simulations require global Poisson-like field solves at 
each time step. This prevents EMGT simulations from scaling as favorably Madduri et al. (2011) 
as VMKT simulations when the number of processing cores is increased at fixed problem size. 
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Likewise, the free energy of perturbations to EMGT equilibria is unknown. Thus, the basie tool 
for studying the stability of EMG equilibria by way of an energy prineiple is unavailable. The pur¬ 
pose of this Chapter is to deseribe a new formulation of eleetromagnetie gyrokineties that enjoys 
properties (I) and (II). The new formulation, whieh we will refer to as the gyrokinetic Vlasov- 
Maxwell (GVM) system, enjoys a local field solve and has an energy principle, while retaining the 
traditional advantages of gyrokinetic theory. 


2.2 The new formulation 


The gyrokinetic Vlasov-Maxwell equations are given by 


dfs 

dt 

IdD 

c dt 
IdB 

c dt 

V D 

V B 


-Lypfs 

(2.1a) 

471 

V X iT- Jgy 

(2.1b) 

-V X E 

(2.1c) 

dTrpgy 

(2.1d) 

0. 

(2.1e) 


fs is the gyrocenter volume form of species s, Vp is the gyrocenter phase space velocity, Lyp 
denotes the Eie derivative along the gyrocenter phase space velocity, Jgy is the gyrocenter current 
density, pgy is the gyrocenter charge density, E, B are the fluctuating electric and magnetic fields, 
and D,H are the auxiliary electric and magnetic fields. The volume form fs is defined by re¬ 
quiring that the number of particles of species s in a region of phase space U be given by fg. 
The gyrocenter phase space velocity is specified by the time-dependent tensor form of Hamilton’s 
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equations, 


i^^gya;f = dKs - CsE ■ dX, (2.2) 

where is the gyrocenter symplectic form, Kg is the gyrocenter kinetic energy, and dX denotes 
the vector line element in the space of gyrocenter positions. The gyrocenter symplectic form is the 
sum of the guiding center symplectic form Cary and Brizard (2009); Burby et al. (2013a) and the 
fluctuating magnetic flux. 


uj^y = ojf - -B ■ dS, 


c 


(2.3) 


where dS is the surface element in the space of gyrocenter positions. The gyrocenter ki netic 
energy is a functional of the fluctuating electric and magnetic fields, and is related to the gyrocenter 
Hamiltonian by = Kg + Cgcj) (an explicit expression for Kg will be given near the end of this 
Chapter). The auxiliary fields D, H are related to E, B by using relations that emerge from the 
Hamiltonian theory developed in Morrison (2013), i.e., the constitutive relations are given by 


D 

H 




B + 4:71 


6B' 


(2.4) 

(2.5) 


where /C(/, E, B) = ^ J fsKg{E, B). 

Following Morrison (2013) the system above constitutes an infinite-dimensional Hamiltonian 
system with dynamical variables /, D, and B, and Hamiltonian functional given by 


'H{f,D,B) = IC{f,E,B) + j P-Ed^X 


( 2 . 6 ) 
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where E = E{f,D,B) is the eleetrie field operator defined implieitly by the equation 


D = E{f, D, B) - 47r^(/, E{f, D, B),B), 


(2.7) 


and P = P(f,D,B) is the gyroeenter polarization operator given by 


P{f,D,B) = —{D-E{f,D,B)). 


( 2 . 8 ) 


The noneanonical Poisson braeket is given by 





- 47re^— ■ dX, d— - 4716^— • dX 
oL) ojs oLf 


fs 


+ Aire 


-• V X -- V X — d X. 

5D 5B 5D 5B J 


(2.9) 


Here is the gyroeenter Poisson tensor, whieh is defined as follows. If z°' is a eoordinate system 
on the gyroeenter phase spaee and a, (3 are 1-forms on the same space, Bf{a, (3) = aal3b{z°'-, 
where {•, is the gyroeenter Poisson bracket. Note that a Poisson bracket for electrostatic gy- 
rokinetics was given in Squire et al. (2013). The complexity of that bracket should be contrasted 
with the relative simplicity of the bracket given here for electromagnetic gyrokinetics. This bracket, 
which has a form akin to that of Morrison (2013), is to our knowledge the first demonstration of 
Hamiltonian structure for any electromagnetic gyrokinetic theory. 


2.3 Origins and comparisons 

We arrived at this electromagnetic gyrokinetic system by modifying the standard variational deriva¬ 
tion of electromagnetic gyrokinetics Sugama (2000); Brizard (2000b,a); Brizard and Hahm (2007); 
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Squire et al. (2013). In the standard approach, a gyrokinetic system Lagrangian is constructed by 
adding a gauge-dependent Xu and Wang (2013) net gyrocenter Lagrangian to a non-relativistic 
limit (known as the Darwin limit) of the free Maxwell field Lagrangian. Applying appropriate vari¬ 
ations to the system Lagrangian then produces the standard equations of EMGT. Roughly speaking, 
adopting a gyrocenter Lagrangian instead of a particle Lagrangian amounts to dropping terms from 
the particle equations of motion. Likewise, adopting the Darwin approximation amounts to drop¬ 
ping terms from Maxwell’s equations. We modified this approach by adding a manifestly gauge 
invariant net gyrocenter Lagrangian Pfirsch and Correa-Restrepo (2004) to the full free Maxwell 
field Lagrangian to produce the system Lagrangian. Thus, in the modified approach, fewer terms 
are dropped from Maxwell’s equations. While dropping these terms as in the standard approach 
would be justified (using the assumption of non-relativistic particles), doing so is not necessary. 
We therefore conclude that the GVM equations are no less accurate than standard EMGT. 

2.4 Computational benefits 

The usual argument for invoking the Darwin approximation in EMGT is that doing so eliminates 
light waves. This may seem to be an especially compelling argument from a computational point 
of view. After all, the presence of traveling waves with phase velocity c leads to a very restrictive 
CEE condition for explicit integration schemes. Therefore, avoiding the Darwin approximation as 
we have done may appear objectionable in a practical sense. 

On the other hand, this numerical argument supporting the Darwin approximation is not as 
strong as it appears. As is evident from the form of the GVM equations given above, avoiding 
the Darwin approximation does not lead to Maxwell’s equations, but Maxwell’s equations in a 
polarized and magnetized medium. Therefore, the light waves supported by these equations do not 
travel at the speed of light in vacuum. 
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It is well known Krommes (1993) that the dieleetrie eonstant resulting from gyroeenter polar¬ 
ization is large, whieh implies that the speed of light is much smaller than c in a gyrokinetic plasma 
(this is consistent with the notion of a so-called “gyrokinetic vacuum”). Using the long-wavelength 
limit of the gyrokinetic dielectric function, as a rough approximation, we find that light 

waves in the GVM equations propagate at the Alfven speed. Thus, the CFL constraint imposed by 
light waves in the GVM equations is not nearly as strict as the usual argument might suggest ^ 

An even stronger case can be made for the computational viability of this new formulation of 
electromagnetic gyrokinetics. We first make the following simple observation. A familiar calcula¬ 
tion shows that if V • T? = dvrpgy and V ■ S = 0 at f = 0, then these equations will also be satisfied 
for all subsequent times. This means that the evolution of the magnetic field and the auxiliary elec¬ 
tric field is completely determined by the Ampere equation and the Faraday equation. Interestingly, 
it can be shown that this property arises as a direct consequence of employing a gauge-invariant 
gyroeenter Lagrangian; the quantity SJ ■ D — dvrpgy is the conserved quantity associated with gauge 
symmetry by Noether’s theorem. 

Now suppose the Ampere and Faraday equations were used to advance D and B in time on a 
computer. Employing a simple explicit scheme, the following steps would have to be taken at each 
time step. (1) Using the constitutive relations, compute E and H from the known values of D and 
B. (2) Compute V x iT and V x E. (3) Using a finite difference approximation for the partial 
time derivative, solve for the new D and B. 

Steps (2) and (3) clearly require only local operations, and so represent nearly embarrassingly 
parallel computations. Again invoking the long wavelength limit, step (1) can also be seen to be 
local. In this limit, there is a simple algebraic relationship between D and E (see Brizard (2013), 
for example) that can be inverted analytically. Thus, the entire field solve step in an explicit time 

* Strictly speaking, it is only light waves that travel perpendicular to the magnetic held that experience a reduced 
propgation speed. Those that travel along the magnetic held lines may still travel near the speed of light in vacuum. 
However, the numerical grids appropriate for gyrokinetic simulations are signihcantly elongated along the held lines, 
which substantially reduces the parallel CFL condition. 




Chapter 2. Hamiltonian formulation of the gyrokinetic Vlasov-Maxwell 

EQUATIONS 


12 


marching scheme for the GVM equations is nearly embarrassingly parallel. Sueh a field solve is 
preferable to the nonlocal Poisson-like solves neeessary in eonventional EMGT, especially when 
performing parallel simulations with very few partieles per proeessing core. 

2.5 Theoretical benefits 

We will now turn from numerieal benefits offered by the GVM equations in order to diseuss their 
analytical benefits. First, we mention the system’s conservative properties. An immediate eon- 
sequence of the GVM Poisson braeket structure is eonservation of the Hamiltonian functional 
(this follows from antisymmetry of the braeket). It is also not diffieult to show that there is 
a eonserved momentum funetional for eaeh rotation or translation symmetry of the baekground 
magnetie field. Finally, there is a large family of eonserved funetionals given by the Poisson 
braeket’s Casimirs. These are functionals C that Poisson commute with every other funetional, i.e. 
VJ^, [C, J^] = 0. Systems of gyrokinetie equations (electromagnetie or electrostatie) with exact 
energy and momentum eonservation laws can also be derived using the standard variational ap- 
proaeh Seott and Smirnov (2010); Sugama (2000); Brizard (2000b,a); Pfirsch and Correa-Restrepo 
(2004); Squire et al. (2013). Indeed, this was the main motivation for developing the standard vari¬ 
ational formulations of gyrokineties. However, variational approaches do not readily produee the 
Casimir invariants (nor has it been shown that the usual variational formulations of FMGT possess 
Poisson braekets and Casimir invariants at all). 

Many of the GVM braeket’s Casimirs are given as follows. Fet 

= Au)f 


( 2 . 10 ) 
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be the Liouville volume form defined by the gyroeenter sympleetic form and introduee the gyro- 
center distribution/wncrion, Fg, where 


fs = Fgng, ( 2 . 11 ) 

then 

Ns ^ 

Ch = y2 hs{Fg)ng ( 2 . 12 ) 

s=l JtQ 

is a Casimir for each function of a single real variable hg. Moreover, any functional of V-Z)—dvrpgy 
is a Casimir, which is one way of seeing that Eq. (2. Id) is satisfied in the Hamiltonian formulation 
of the GVM equations. 

Another advantage the Poisson bracket formulation of the GVM equations provides, which a 
variational formulation does not, is immediate access to the theory of dynamically accessible vari¬ 
ations Morrison and Pfirsch (1989) (see also Morrison (1998); Andreussi et al. (2013)). Suppose 
we perturb a GVM equilibrium by switching on a small time-dependent term in the Hamiltonain, 
i.e. Ti ^ y. + where ST-Lt is a time-dependent functional that is non-zero only in a brief 
interval of time after f = 0. Using the Poisson bracket, we can give an energy principle for as¬ 
sessing the stability of this perturbation in the limit where the kick caused by switching on SHt is 
infinitesimal. 

In this limit, and accounting for the fact that the perturbation is generated by altering the Hamil¬ 
tonian, we find that the perturbed distribution function, auxiliary electric field, and magnetic field 
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must have the form 


6fs = -L^Js 

6D = -47r J(^, /) + 47rc V X /3 


(2.13) 


(2.14) 


6B = —47rcV X a, 


(2.15) 


where a, /3 are arbitrary vector fields on configuration space, the phase space fluid displacement 
vector is determined by Hamilton’s equations, 


i^^r = dXs + 4:7resOL ■ dX, 


(2.16) 


with Xs an arbitrary function on gyrocenter phase space, and f ) is the gyrocenter current den¬ 

sity generated by fiducial gyrocenters with phase space velocity and distribution fg. Appealing 
to the general theory of dynamically accessible variations (see e.g. Morrison (1998)), our pertur¬ 
bation will be stable if the free energy functional 5‘^F{o!., j3, x) is positive whenever 5fs, 5D, and 
5B dse not each zero. The free energy functional is defined by 




(2.17) 


where the functional S = Yl,g f Xsfs + f ct -D dfX + fj3-B dfX. Physically, 5‘^F is the second- 
order change in the energy functional "H produced by our perturbation. In fact, 5‘^F functions as 
the (conserved) Hamiltonian of the linearized GVM equations. 
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We find that 5‘^F can be written in the form 

= j Sf. 

+ — [ (sD-dE 

SttJ \ 


+ 6B-6H] FX. 


(2.18) 


Here X in a subscript denotes the X-component of a velocity field on phase space. The variations 
5Ks, 6E, and 6H are given by 



(2.19) 

6E = e-^ [6 D]+r][SB] 

(2.20) 

6H = 7]^[6D]+ fi-^[SB]. 

(2.21) 

where the linear operators e, /i, and r] are given by (cf. Monison (2013)) 


£ = 1 — AtT - 

6E6E 

(2.22) 

^ , 5^ 1 C ,2 52/c S^JC 

^ ~ 5E6B^ 5B6E 

(2.23) 

SBSE- 

(2.24) 


In principle, an energy principle for electrostatic gyrokinetics analogous to this one could be de¬ 
rived using the Poisson bracket given in Squire et al. (2013). However, the authors of that Refer¬ 
ence deemed the electrostatic gyrokinetic Poisson bracket too complicated to be practically useful, 
and so did not attempt deriving an expression for 5‘^F. 

We have used this expression for S‘^F to prove that, in the long wavelength limit, the thermal 
equilibrium state in a uniform background magnetic field is stable. In this case, the gyrocenter 
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kinetic energy is given by 


K 


1 2 T 1 2 

-mvu + uJcJ - me 

2 II 2 


\c Bo 


E xb 



2 




(2.25) 


where J is the gyroaction, ujc is the signed gyrofrequency, = B — bb ■ B, and Bo is the mag¬ 
nitude of the background magnetic field. This expression agrees with that given by Krommes in 
Krommes (2013) in the absence of magnetic fluctuations. The linear response functions r] 

are therefore given by the constant matrices 


47rc 

1 + 2 (1 bb) 

(2.26) 

n 


1 — 47r/9(l — bb) 

(2.27) 

0, 

(2.28) 


where /? = E. is the plasma {3 and {■)s denotes the velocity space average. Using these 

expressions and the assumption of thermal equilibrium, a straightforward, but tedious calculation 
leads to the following form for 6‘^F, 




Stt J 

[ 5D-e-HDd^X 


Stt j 

[ 6B^ ■ /i-^ • SB^d^X 


Stt j 

[{1- ^'KnTjBDbBl d^X, 

(2.29) 


where n = is total gyrocenter number density. As long as dvr/? and AimT / B^ are each 

less than 1, a condition that is generally satisfied, 5‘^F is manifestly non-negative, which implies 
linear stability. 
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2.6 Concluding remarks 

The Hamiltonian formulation of the GVM system given in this Chapter is eompletely determined 
by two key quantities, the gyrocenter kinetie energy Kg, and the guiding eenter sympleetie form 
L;f. Suppressing speeies labels, the gyroeenter kinetic energy is given explicitly to second order 
in the amplitude of the fluctuating fields, es, by 


K{E, B) = - es (£) + , d (£)) 



(2.30) 


where R is the infinitesimal generator of gyrophase rotations times the local gyrofrequency, I is the 
inverse of the Lie derivative Lr, angle brackets denote gyroangle averaging, and Q = Q — (Q). In 
standard guiding center coordinates, Lr = Uc^, where 9 is the gyrophase, which means I amounts 
to an antiderivative in gyrophase. It can be shown that the second-order gyrocenter kinetic energy 
has the same general form as Eq. (129) in Brizard and Hahm (2007). The relevant correspondences 
between our symbols and those of Brizard and Hahm (2007) are £ —Ki, -H- 55 -H- 

AT, and ^ Lii[Ti + Ti). 

From this expression, it is clear that the gyrocenter kinetic energy is determined by the three 
quantities i, and 55. denotes the guiding center Hamiltonian truncated at some desired 
order in p/L. The function i and the 1-form 55 are defined in terms of any choice of the guiding 
center Lie generators as follows. Decompose the guiding center transformation Toc ■ TQ ^ TQ 
as Tgc = r 2 o n, where 


Ti = exp(Gi) 


(2.31) 


T2 = ■■■o exp(G3) o exp(G2) = exp(G2), 


(2.32) 
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and the Gk are the guiding eenter Lie generators. The leading-order guiding eenter transformation, 
Ti, must be handled carefully in gyrokinetics because the fluctuating fields are allowed to have 
short perpendicular wave lengths. The 1-form 

+ iG,U{La,))B ■ dS, (2.33) 

where the function U{x) = e“'^/^sinh(a:/2)/(a;/2), represents the perturbation to the guiding 
center Lagrange 1-form produced by the fluctuating electromagnetic fields. The function 

5n = e{exp{-L^JiG,U{LG,) + iG,U{LG,))E ■ dX (2.34) 

represents the perturbation to the guiding center kinetic energy caused by the same fields. The 
function 


f = 5S(l/s5') - 5 H , (2.35) 

where is the unperturbred gyrocenter phase space velocity. 

The Hamiltonian structure of the GVM equations reproduces that of the Vlasov-Maxwell sys¬ 
tem Morrison (1980, 1982); Marsden and Weinstein (1982) under the substitutions 

K —)■ (2.36) 

—)■ m dx* A dvi- (2.37) 

It is also interesting to compare [•, •] to the bracket given in Morrison (2013). The only significant 
difference comes from the manner in which the inductive electric field is built into the ki netic 
equation. 
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Finally, we note two possible direetions for future researeh. (1) It may be useful to identify 
a Poisson braeket for electromagnetic gyrokinetics in the Darwin approximation, i.e. standard 
EMGT. The gyrokinetic Vlasov-Darwin equations are somtimes also referred to as the gyroki¬ 
netic Vlasov-Poisson-Ampere equations Sugama (2000). A Hamiltonian formulation of the non- 
gyrokinetic Vlasov-Darwin equations has already been given in Krause et al. (2007). (2) It seems 
likely that the bracket and Hamiltonian given in this Chapter will provide the Hamiltonian structure 
for the oscillation center Vlasov-Maxwell equations with appropriate substitutions for K and 
If this were true, then the benefits that our bracket brings to electromagnetic gyrokinetics could be 
extended to certain kinds of laser-plasma interactions. 
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Chapter 3 

Boundary terms and Poisson brackets 

3.1 Introduction 

The purpose of this Chapter is twofold. The first is to deseribe a systematic procedure to pass from 
an Euler-Poincare formulation Holm et al. (1998) of a physical system with advected parameters 
and dynamical fields to a Poisson bracket formulation for that same system. The second is to apply 
this procedure to the Euler-Poincare formulation of several reduced plasma models, including the 
gyrokinetic Vlasov-Maxwell system presented in the previous Chapter. The value of a general 
Eluer-Poincare—> Poisson procedure stems from the fact that implementing approximations within 
the Eagrangian formalism is a well-developed art, whereas finding approximations that respect the 
Jacobi identity is much more subtle. 

To pass from an Euler-Poincare formulation with a regular Eagrangian to a Poisson bracket 
formulation, all that is necessary is the procedure outlined in Holm’s paper on Euler-Poincare 
theory Holm et al. (1998), which consists of two steps. Eirst one passes from the parameterized 
Eagrangian description to a parameterized Hamiltonian description using the Eegendre transform. 
Then one applies the theory developed by Marsden in his paper Marsden et al. ( 1 984) on the Hamil- 
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Ionian side of semi-direct product theory to “give life” to the advected parameter and identify a 
Poisson bracket for the total system. 

When the Lagrangian is degenerate, the identification of an appropriate Poisson bracket is 
more difficult; the conventional Legendre transform-based technique fails. Nevertheless, using 
the theory of Dirac constraints, CendraCendra et al. (1998) and Squire Squire et al. (2013) have 
derived brackets for specific Euler-Poincare systems with degenerate Lagrangians. Both Cendra 
and Squire have rederived the Morrison-Marsden-Weinstein bracket. Squire has derived a bracket 
for the Gyrokinetic Vlasov-Poisson equation. In principle, the methods used by these authors 
could be adapted to suit our needs, but this is not the path we will follow. Instead we will develop 
a technique that compliments those of Cendra and Squire. The motivation behind “reinventing 
the wheel” is our desire to simplify the Cendra-Squire approach in such a way that clarifies how 
advected parameters become dynamical variables and eliminates the need to explicitly introduce 
the Dirac theory of constraints. 

The technique we will develop is insensitive to the degeneracy of the parameter-dependent 
Lagrangian. In this sense it is similar to the method of Cendra and Squire. However, we will not 
employ the Dirac theory of constraints. Instead, we will apply a technique rooted in a careful 
analysis of the boundary terms that appear when varying an action functional without keeping 
endpoints fixed. In detail, our method consists of the following steps. 

step 1 — Identify an Euler-Poincare formulation for the system under consideration. In par¬ 
ticular, identify a parameter-dependent Lagrangian, La TQ x TG ^ M., where Q is the space 
of dynamical fields, G is a Lie group (usually a diffeomorphism group), and the parameter a is an 
element of a vector space V* upon which G acts. 

step 2 — Eliminate the parameters by introducing a Lagrange multiplier. This method is de¬ 
scribed in Cendra’s Lagrangian reduction by Cendra et al. (2001). The result of this simple 
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step will be a Lagrangian C : TQ x TG x T{V xC*)—that embeds the original dynamies in 
a slightly larger spaee. 

step 3 — Identify the submanifold Vo C TQ x TG x T{V x V*) that serves as the aug¬ 
mented system’s phase spaee. Note that this step involves analyzing the initial value problem 
assoeiated with £’s Euler-Lagrange equations. Possible gauge symmetries and degeneraeies of 
the Lagrangian make this step non-trivial in general. 

step 4 — Identify a Poisson braeket h ■]-Po and Hamiltonian TL-Po on the augmented system’s 
phase spaee using the boundary sympleetie form methodology deseribed in Marsden et al. (1998). 
This braeket is essentially a Pierles braeket. 

step 5 — Observe that the Pierles braeket and Hamiltonian on the augmented phase spaee are 
invariant under the aetion of the semidireet produet S = G x V. Perform Poisson reduetion using 
this symmetry, thereby identifying the Poisson braeket on the redueed phase spaee V = VojS. In 
this step, the Lagrange multiplier will be eliminated by the Poisson reduetion, meaning V ean be 
thought of as the physical phase space for the system under consideration. 

In what follows, we will illustrate this technique by applying it to a number of examples. 
These include (i) the (generalized) Vlasov-Poisson system, (ii) the Vlasov-Darwin system, (iii) 
the gyrokinetic Vlasov-Maxwell system, (iv) the neutral Vlasov model Tronci and Camporeale 
(2015), and (v) force-free electrodynamics Gralla and Jacobsen (2014). 
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3.2 The generalized Vlasov-Poisson system 

3.2.1 Step 1: Euler-Poincare formulation 

Let P be a 2iV-dimensional symplectic manifold with symplectic form u = —note that P is 
not necessarily a cotangent bundle and u is not necessarily a canonical symplectic form. Let V* = 
^ 2 Af(p) denote the space of 2iV-forms on P and the group of diffeomorphisms 

of P, respectively. A typical element of G will be denoted g G G while a typical element of V* 
will be denoted f E V*. 

The generalized Vlasov-Poisson system is defined by the parameter-dependent Lagrangian, 
Lf^ : TG -E R, given by 

Lfoi9:9) = J i9Jo)^{g o 9~^) -'H{9*fo), (3.1) 

where "H : V* —)■ R is the generalized Hamiltonian functional. The generalized Vlasov-Poisson 
dynamics follow from this Lagrangian by applying Hamilton’s principle to the following action 
functional. Let ^(G) and g G ^(G) denote the space of paths in G and a typical path in G, 
respectively. The action functional Sf^ : ^(G) —)■ R is given by 

Sfois)= [ (3.2) 

Jti 

The Euler-Lagrange equations associated with this action functional can be derived as follows. 
Let 0 = X(P) and G 0 denote the space of vector fields on P and a typical vector field, respec¬ 
tively. Notice that 


Lfoig^g) = ^{g°g \g*fo), 


(3.3) 
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where £ : g x C* ^ R is given by 


mj) = jjm-nf)- 


(3.4) 


It follows that the first (fixed-endpoint) variation of the aetion funetional is given by 




fp ^ ^ “ I7 dt 

f f dSi ^ Si Si \ ^ ^ , 


rt2 


Iti JP 


^ ® fit) + L^(t) (i? ® fit)) - d i)(|(t)) - 


SH 

w 


CL 


fit) (^i^p)a; - ■ r/(t) dt, 


fit) ) ■ riit) dt 
(3.5) 


where 


^it) = g(t) og(t)"^ 

(3.6) 

rfit) = Sgit)og{t)-^ 

(3.7) 

fit) = sit) Jo, 

(3.8) 


and the functional derivatives are evaluated at (^(t), fit)). The Euler-Poincare equations are there¬ 
fore 


k(t)^ 



which should be augmented with the equation 


fit) = git)*fo. 


(3.9) 


(3.10) 
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Note that the last equation implies f{t) = 

3.2.2 Step 2: Introduction of the Lagrange multiplier 

Note that the equation f{t) = from the previous section does not follow from the Euler- 

Lagrange equations associated with the parameter-dependent Lagrangian L because fo is treated 
as merely a parameter in Hamilton’s principle. We can formally remedy this issue as follows. 

Let V = V** and x ^ ^ be the dual to V* and a typical element of V, respectively. The 
space V is properly the space of continuous linear functionals on V*, which is naturally the set 
of distributional functions on P. Define the parameter-independent Lagrangian C : TG x T{V x 
V*) —)■ R by 


^9,9,xJo,xJo) = LfXg.g) 



(3.11) 


When Hamilton’s principle is applied to the augmented action functional, S : ^(Gx Cx C*) —)■ R, 
given by 

S{s,X,fo)=[ C{g{t),g{t),x{t),fo{t),x{t),fo{t))dt, (3.12) 

Jti 

the resulting Euler-Lagrange equations are given by 


fo{t) = 0 

x{t) = g{tT (^i?(^(t)) - 

= dyr, 


(3.13) 

(3.14) 

(3.15) 


where the functional derivatives are evaluated at g(t)*/o(f)- We have thus succeeded in embedding 
the Euler-Poincare equations along with the advection equation into a larger system. We will refer 
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to the system defined by Eqs.(3.13), (3.14), and (3.15) as the augmented generalized Vlasov- 
Poisson equations (AGVP equations, for short). 

3.2.3 Step 3: identification of the augmented phase space 

We will now study the initial value problem assoeiated with the augmented generalized Vlasov- 
Poisson equations. In partieular, we would like to identify a submanifold, Vo C TG x T{y x V*), 
sueh that the AGVP equations define a first-order (infinite-dimensional) ODE on Vo- In order to 
aeeomplish this task, we will merely rearrange the AGVP equations given in the previous seetion 
into the form of a first-order equation, and then deduee the allowed set of initial data. 

As they were written in the previous seetion, the AGPV equations are nearly expressed as a 
first-order system. In order to aehieve the desired form, we re-write Eq. (3.15) in terms of g(f) and 
substitute Eq. (3.15) into Eq. (3.15), giving 

fo { t ) = 0 (3.16) 

W) = gM* - ^) (3-17) 

g { t ) = Xmog { t ), (3.18) 

where the funetional derivatives are evaluated at g(f)*/o(f), whieh is elearly a first-order system 
of equations in the variables (/o, x, g). That is, there is a veetor field Y on G x V x V* sueh that 

(gW, g(f), XW, X(f), /o(f)) = y (g(f), X(f), fo { t )). (3.19) 

The triple {foiXid) belongs to the set G x V x V*, whieh can naturally be identified with 
the graph of the vector field Y on G x V x V* that is defined by Eqs. (3.16), (3.17), and (3.18). 
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Therefore the AGVP equations naturally define a first-order system of ODEs on the submanifold 


Vo = {{g, 9 . X, fo, X, fo) eTGx T{V x C*) I {g, g, x, fo, X, fo) = Y{g, x, fo)} ^ G x V x V*. 


(3.20) 


Interestingly, Eq. (3.18) implies that there is an invariant subset of Vo given by 

Vo = {{g,g,x,fo,xJo) e Vo\g e DifL(P)} ^ DifL(P) x C x C*, (3.21) 

where DifCj(P) is the set of sympleetie diffeomorphisms of P. However, Vo is not in one-to-one 
eorrespondenee with all solutions of the AGVP equations, whereas Vo is. Therefore we will regard 
Vo as the phase space for the AGVP equations. 

3.2.4 Step 4: derivation of the boundary sympleetie form on Vo 

Because Vo is a valid phase space for the AGVP equations, the AGVP dynamics formally define a 
time-independent flow map '^t ■ Vo = G x V x V* ^ Vo = G x V x V*, which is characterized 
by the relations 


(3.22) 

j^U9,xJo) = YiU9,xJo)), (3.23) 

where Y is the vector field on G x V x V* defined by Eqs.(3.16), (3.17), and (3.18). We can 
therefore define a mapping Sol : Vo ^ ^{G x V x V*) given by 


Sol(^,x,/o)(f) = dt-h{g,x,fo)- 


(3.24) 
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The mapping Sol can be used to pull back the augmented action functional S to the augmented 
phase space Vo, thereby defining the restricted augmented action 

= SorS. (3.25) 

By examining the exterior derivative of the restricted augmented action, we can identify a sym- 
plectic form, and therefore Poisson brackets, on the augmented phase space Vo- This can be seen 
using the following formal manipulation. 

Let {g, X, fo) £ Vo be an arbitrary point in the augmented phase space. Let 

fo{t)) =dt-ti{9,X,fo)- (3.26) 

Because the AGVP dynamical equations imply fo{t) = fo, the restricted augmented action evalu¬ 
ated at {g, X, fo) is given by 

Spo{9,X,fo) = [ Lf^{g{t),g{t))dt= [ g{t)Jo), (3.27) 

Jti Jti 

where $,{t) = g{t) o g(t)~^. Using the AGVP equations of motion, we can therefore write 

dSp^ = - S, (3.28) 

where S is a one-form on Vo given by 

^{9, X, fo) [59, 5x, 5fo] = j x5fo + ^{^9 ° 9~^) 9*fo- 


(3.29) 
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Differentiating Eq. (3.28) in t 2 (and remembering that depends on ^ 2 ), we obtain 


dSr^ = LyE ^ iydE = -d (h(F) - , 

where the funetional : Vo ^ ^ is given by 


(3.30) 


^Void, X,fo) — ^ {^sn/5f{g,fo),9*fo) ■ (3.31) 

Equation (3.30) immediately implies that the 2-form oo-p^ = —dS is preserved by the AGVP flow, 
Moreover, beeause it is not hard to show that u-p^ is non-degenerate, the AGVP equations ean 
be written in Poisson braeket form, i.e. given a functional V : Vo ^ 


V = [V, Hp^] 


Vo-, 


(3.32) 


where 


Hp^{9,xJo) = ^{Y)-Sp^ 

= nig Jo) (3.33) 

is the augmented system’s energy functional and [•, - is the Poisson bracket obtained by inverting 
the two-form ojp^. 

We will conclude this section by deriving an explicit expression for the bracket [•) ^Vo- Eirst 
observe that if ; Po —^ is a functional, then the associated Hamiltonian vector field Yp, i.e. 
the vector field on Vo that satisfies 


iy^uJVo = 


(3.34) 
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is related to the Poisson braeket aeeording to the formula 

Lx^g = [g,J^]p^. (3.35) 

Beeause the Lie derivative Lxj,g would be easy to ealeulate if Xjr were known, we will eompute a 
general expression for Xjr- To this end, it is useful to observe that any integral curve of the vector 
field Xjr is a critical point of a phase space variational principle. That is, an integral curve of Yjr is 
automatically a critical point of action functional Ax : x L x L*) —)■ M given by 

^A9,xJo) = [A9it)Xt)Jo{t))[{At),At)Jo{t))] - X{g{t),x{t)Jo{t))^ dt. (3.36) 

Therefore we can derive an expression for Yx by varying the action given in Eq. (3.36). For the 
sake of varying T w.r.t. g, we introduce the convention that if if is a functional on Vo, then 
5H/5g{g, x, fo) is the unique 1-form on P that satisfies 

dH{g,xJo)[5g,0,0] = j (^^ig,X,fo)®g*f<^-i5gog-^). (3.37) 

The first (fixed-endpoint) variation of Ax is given by 


SAxig, X, fo)[5g, 5x, 5fo] = 



5x + 


Y - x(i)) mt) 
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where 


^{t) = g{t) O g(t) 1 (3.39) 

r}{t) = 6g{t) o g(t)"\ (3.40) 

and funetional derivatives are evaluated at (g(t), x(t), /o(^))- It follows that Yjr is speeified by the 
relations 


g{t) = uj ^ (^SJ^/Sg + ° S(^) (3-41) 

Xit) = gity (^6J^/6g + ^ ^ - 5T/5fo (3.42) 

Ut) = ST/^X, (3.43) 


where oj ^ denotes the inverse of the linear map X i^ce. This formula for Yjr proves that the 
two-form —dS is non-degenerate. Moreover, we can now write down the Poisson bracket [t -jPo 
using Eq. (3.35). The result is 


[J^,G]vo 



(b7_ ^ 



9*fo + 


SfoSx 



Sfo 




(3.44) 


where B denotes the Poisson tensor assocaited with the symplectic form u, i.e. given 1-forms on 
P, a and (3, 


B{a,l3) = u(u \a),u ^(/ 3 )). 


( 3 . 45 ) 
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3.2.5 Step 5: perform Poisson reduction to obtain bracket on physical phase 
space 

We have now identified the Hamiltonian, 


'Hr^{9,xJo)=ngJo), ( 3 . 46 ) 

and Poisson bracket for the augmented generalized Vlasov-Poisson system. The Jacobi identity is 
satisfied because the bracket has been obtained by inverting the symplectic form —dS. The modi¬ 
fier “augmented” is appropriate because the dynamical variable x has no direct physical meaning. 
On the other hand, the variables g and fo together comprise an element of the physical Lagrangian 
(as opposed to Eulerian) phase space; g gives the configuration of particles in the single-particle 
phase space P and fo gives the reference phase space density. The purpose of the additional vari¬ 
able X is to extend the Lagrangian phase space just enough to allow for a non-degenerate Poisson 
bracket. 

The appearance of the variable x perhaps seems awkward at this stage. However, observe 
the following. The set C is a Lie group under addition that is a symmetry group for the AG VP 
equations. Specifically, for each (5x G V, we can define a mapping ■ Po ^ Po given by 


Tsxia, X, fo) = {g,x + ^x, fo), 

which clearly satisfies the defining properties of a group action, 

^5x1+15x2 ~ '^Sxi ° ^<5x2 


( 3 . 47 ) 


( 3 . 48 ) 


To — idp^. 


(3.49) 
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This group action leaves the augmented Hamiltonian and the augmented Poisson bracket 
[•, ■]'p^ invariant in the sense that 


= Hv, (3.50) 

(3.5i) 

for arbitrary 5x ^ V and funetionals Q. These properties are quiek to verify. Therefore we 
ean define a Hamiltonian and Poisson braeket on the first reduced phase spaee Vl = Vo/V, whieh 
is naturally the physieal Lagrangian phase spaee, G xV*. The Lagrangian Hamiltonian, 'H'p^ : 
Vl —)■ R, is simply given by 


nv,{gJo)=n{gJo). (3.52) 

The Lagrangian Poisson braeket, [•, is defined in terms of the projeetion map ttl ■ Vo ^ Vl 
given by 


7^L{9,x,fo) = {gjo)- 


(3.53) 


We have 


7rl[d,e]p, = [7rld,7rl&]p^, 


(3.54) 
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which defines [•, unambiguously beeause Ts^ leaves the Poisson braeket invariant. In partieu 


lar, 


[5, <&]vl{9Jo) = 


(3.55) 


L 


S(5 



(3.56) 


where x is arbitrary and the funetional derivatives are evaluated at {g, fo). It does not matter whieh 
X is ehosen beeause 



= [Ti^TTld, Ti^7rl(5]rA9, X, fo) (by Eq. (3.193)) 

= [TTld, 7r*L^]pA9, X, fo) (by Eq. (3.53)). (3.57) 


The relation given in Eq. (3.54) shows that the Eagrangian braeket automatieally satisfies the Jaeobi 
identity. Thus, the awkwardness introdueed by the additional variable x is only apparent; we have 
obtained a physical Hamiltonian formulation for the generalized Vlasov-Poisson in Eagrangian 
labeling by reeognizing that the augmented Hamiltonian and braeket are independent of x- This 
Hamiltonian formulation is “physieal” in the sense that all dynamieal variables are physieally 
signifieant. 

The set G is a symmetry group of the generalized Vlasov-Poisson system in Eagrangian label¬ 
ing. Speeifieally, for eaeh h e G, we ean define a mapping Rh : Vl ^ Vl given by 


Rh{9,fo) = {9°h, h*fo), 


(3.58) 
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that satisfies the properties of a right group aetion, namely 

Rhioh2 = Rh2 ° Rhi (3.59) 

Ridp = idpi- (3.60) 

The Lagrangian Hamiltonian and Poisson braeket are eaeh invariant under this group aetion, in the 
sense that 


RlRvp = Rvp (3.61) 

Ki[r,gw = \Kr,RlQW (3.62) 

The invarianee of the Lagrangian Hamiltonian is quiek to verify. The invarianee of the Poisson 
bracket follows from the identity 

(g, /„) = ^(9 ° fc. fc*/.). (3.63) 

Therefore, the Lagrangian Hamiltonian and Poisson bracket define corresponding quantities on the 
Eulerian phase space, Ve = Vl/G ’^V*. The Eulerian Hamiltonian is given by 

'HvM)=nf)- (3.64) 

The Eulerian Poisson bracket is defined in terms of the Eulerian projection map, tte ■ Vl ^ Re, 
given by 


'^e{,9-i fo) 9*fo) 


(3.65) 
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which satisfies the important property tte o Rh = We have for functionals F, G : —)■ M, 


[F,G]p,(/) 


f B 


Vi/' ^ 


6F 6G 


\ 


/, 


(3.66) 


where g and fo are any group element and 6-form that satisfy / = g^^fo, and is the Poisson 
bracket associated with the symplectic form u. Just as earlier, this bracket automatically satisfies 
the Jacobi identity. However, while this fact for the Lagrangian bracket is perhaps not well-known, 
here it comes as no surprise; [•, - is none other than a Lie-Poisson bracket. 


3.3 The Vlasov-Darwin system 

3.3.1 Step 1: Euler-Poincare formulation 

We will use the following notation. 

Particle configuration space — Let Q = be the single-particle configuration space with 
metric tensor (•, •) and associated hodge star *. Typical elements of Q will be denoted q G Q. The 
codifferential on Q will be denoted S, the Laplace-deRham operator will be denoted A = (5d -f d(5, 
and the Laplace-deRham Green operator will be denoted G. We will also make use of the trans¬ 
verse and longitudinal projection operators 11^ = (5Gd and Ul = dG6. 

Particle phase space — The set TQ will serve as the single-particle velocity phase space. A 
typical element of TQ will be denoted Vq G TgQ. The map tt : TQ ^ Q will denote the tangent 
bundle projection. Let F ; TQ —)■ T*Q be the diffeomorphism given by Vq i—)■ {vg, ■). The symbol 
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will stand for F*6'o, where 6o is the eanonieal 1-form on T*Q. 

Function spaces — The sets Diff(TQ), C°°{TQ), and den{TQ) denote the diffeomorphism 
group of TQ, the spaee of smooth funetions on TQ, and the the spaee of distributional densities 
on TQ. Typieal elements of these spaees will be denoted g G Diff(TQ), x G C°°{TQ), and 
/ G den(TQ). The sets for integer k > 0 are the fc-forms on Q. For our purposes, the 

0-forms, 1-forms, and iV-forms are the most important. Typieal elements of the latter will be 
denoted 0 G A G and A G The sets X{TQ) and X{TQ)* are the veetor 

fields and 1-form densities on TQ. If © is any spaee, ^(©) will denote the spaee of paths in & 
parameterized by the time interval [Q, ^ 2 ]- If the symbol s is used to denote a typieal element of 
6 , we will use a bold version of the same symbol to denote a typieal path in 6, i.e. s G ^(©). 

Multi-species objects —Let Ns be the number of plasma speeies. Set G = Diff(TQ)^% 
V = , and V* = den(TQ)^G We will denote typieal elements of these spaees with 



A typieal element of g will be denoted ^ G g. 

SugamaSugama et al. (2013) gives the following Euler-Poincare formulation of the Vlasov- 
Darwin system (which he calls the Vlasov-Poisson-Ampere system). Define the parameter- 
dependent Lagrangian, Tf : TG X T(D°(Q) x D^(Q) x D^(Q)) —)■ M, given by 

Jo 



rris^ + ids o S'/) - gs*fo,s {Ks + es7i*(f)) 



(3.67) 




Chapter 3. Boundary terms and Poisson brackets 


38 


where Ksivq) = rus {vq,Vq) /2. The Vlasov-Darwin system of equations in Lagrangian labeling 
then follow from Hamilton’s prineiple applied to the aetion funetional x x 

X —)■ R given by 

0, A, A) = f Lj^{g{t),g{t),(l){t),A{t),X{t),<p{t),A{t),X{t))dt. (3.68) 

J 

The Euler-Lagrange equations assoeiated with Sugama’s Lagrangian are given by 




(3.69) 


4 1 

MA(t) = — ^ * u{is{t), fs{t)) -f -d * X{t) 

C “ c 

S = 1 

(3.70) 


Ns 



Sd(j){t) = dTT^e^ * n{fs{t)) 

S = 1 

(3.71) 


o' 

(3.72) 

where (t) = gs{t) 

° Ssit) \ fsit) = gsit)Jo,s, the quantities 



Os = 'fTT's'd + —n*A{t) 
c 

(3.73) 


Hs = Ks + es7r*(p{t), 

(3.74) 

and the operators u 

: X{TQ) X den(Tg) ^ n\Q) and n : den(Tg) ^ n^Q) 

are given by the 

fiber integrals 

/)(?)=/ kf 

Jw l((j) 

(3.75) 


n{f){q) = [ /• 

Jtt l(q) 

(3.76) 


We will refer to u{^, f) as the partiele flux 2-form and n(/) as the partiele spatial density 3-form. 
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We will not use Sugama’s variational formulation as our initial Euler-Poineare formulation for 
the Vlasov-Darwin system. Instead we will insert the elliptie equations for the potentials baek into 
Sugama’s Lagrangian, thereby obtaining a new Lagrangian : TG —)■ M. A straightforward 
ealeulation shows that L has the simple expression 


Lf^{9,9) = i{9°9 ^,9*fo), 


(3.77) 


where 


«(«. /) = (E / ■!>) K*>) - 


(3.78) 


s=l 


Here we have introdueed the Hamiltonian funetional "H ; C* —)■ 




(3.79) 


the potential operators A : g x V* ^ and $ : C* —)■ fl°(Q), 


A(lf) = —GnTJ(iJ) 

c 

$(/) = 47rGp(/), 


(3.80) 

(3.81) 


and the eharge and eurrent density operators p : C* —)■ and J : g x V* ^ 


Ns 


P(f) = *n(/,) 


(3.82) 


5=1 

Ns 


f) = ^es*u{^sjs)- 

5=1 


(3.83) 
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We will now verify direetly that our parameter-dependent Lagrangian reproduees the 
Vlasov-Darwin equations. Because the equations for the potentials are satisfied by construction, 
we will merely verify that the Euler-Poincare equation associated with reproduces Eq. (3.69). 
The Euler-Poincare equation follows from Hamilton’s principle applied to the action functional 
: ^(G) —)■ M given by 

‘5/„(g) = / Lj^{g{t),g{t))dt. (3.84) 

Jti 

Varying this action, we obtain the general Euler-Poincare equation given originally by Holm, 


d6i 5i 


^Js 


(3.85) 


In order to calculate the functional derivatives appearing in this expression, we will first express 
the reduced Eagrangian i in the form 

((I /) = E / + T h. J(l /)) - «(/). (3-86) 

s=l dTQ 

where (•, •) denotes the natural integration pairing of differential forms 


{a, (3) = / a A 
Jq 


(3.87) 


Next we introduce the linear operators : g —)■ f2^((5) and given by 


(3.88) 


= Ai f) 


(3.89) 
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The adjoint operators (J|)j : —)■ C°°{TQ) and {J^)l : —)■ X{TQ)* defined by the 

relations 


Ns 


a 


a 


S=1 JT'Q 

s=l ''T'Q 


(3.90) 

(3.91) 


are readily found to be given by the formulae 


= es7r*a{Q 
{J-^)l{a) = es7r*a O /« 


(3.92) 

(3.93) 


Finally, we eompute the Freehet derivative of i. 


"p 0 /.) . 5i, + («,) - ^ j i/ + i (.4(1, /), J’m) + J (.4(1, /), Jj(i/) 

= E (“•’’ + i(.^^)I(.4({, /))) . Si. + + i( Jf)!(.4(e /)) - SS 

= y ® /» + ((ir‘.4(|, /) 0 • |55, + yn,i?(5,) + ^w'Ad, f)(i.) - Sf, 

(3.94) 


from whieh the funetional derivatives ean be quiekly extraeted, giving 

5i e ~ ~ 

— = rris^ ® fs + —7r*^(^, /) 0 fs (3.95) 

O^g c 

^ = rUgi^iQ + —n*A{i, f){Q - 
djs c df 


(3.96) 
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Note that in computing the Frechet derivative, we have made use of the fact that the operator GUt 
is self-adjoint. Equation (3.85) therefore reduces to 

/(f))) + ^ (m^i^ -f /(f))) + = 0, (3.97) 

which is readily verified to be equivalent to Eq. (3.69). 


3.3.2 Step 2: Introduction of the Lagrange multiplier 


We now define the parameter-independent Lagrangian C : TG x T(y x V*) —)-Mby 


Ns . 

^9,9,xJo,X,fo) = Lj^{g,g) + '^ Xsfo,s- (3.98) 

s=i Jtq 


When Hamilton’s principle is applied to the augmented action functional, S : ^(G x Cx 1/*) —)■ M, 
given by 


‘5(g,X,/o) = / C{g{t),k{t),x{t)Jo{t),ic{t),fo{t))dt, (3.99) 

Jti 


the resulting Euler-Lagrange equations are given by 


= 0 

Xs{t) = = gsity (^Osi^sit)) - 



(3.100) 

(3.101) 

(3.102) 


where 


Os = rusd + -^'n*A{i{f), /(f)), 


(3.103) 
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and the functional derivatives are evaluated at g(t)*/o(^)- The parameter-independent Lagrangian 
C therefore succeeds at embedding the Vlasov-Darwin dynamics, including those of the distribu¬ 
tion function, into a larger system. We will refer to the system defined by Eqs. (3.100), (3.101), 
and (3.102) as the augmented Vlasov-Darwin equations (AVD equations, for short). 

3.3.3 Step 3: Identification of the augmented phase space 

In order to identify a phase space for the AVD equations, we must identify a submanifold of TG 
upon which the Euler-Lagrange equations associated with the augmented Eagrangian C define a 
system of equations that are first order in time. As the Euler-Eagrange equations are written in the 
previous section, the AVD equations are not written as a first order system, and so there is some 
work to do. 

The most non-trivial AVD Euler-Eagrange equation is the one that takes the form of a time- 
dependent Hamilton equation: 

(ras'd + ^7r*A{i{t), f{t))^ + ^ (nis^ + = 0. (3.104) 

This equation appears to implicitly relate the Eulerian velocity field ^ (f) with its time derivative 
^{t) and the distribution function f{t); the time derivative of ^{t) appears as a result of the identity 

4.4(|((),/(()) =.4(|((),/(()) -.4(|((),%,/(()). (3.105) 

Eortunately, this implicit relationship is not quite as complicated as it seems. It turns out that 
Eq. (3.104) implies ^s(t) must be a second-order vector field for each t. Therefore ^(^(f), /(t)) 
and u{^s{t)-i fsil)) can be expressed in terms of the free-streaming vector field Xo, i.e. the unique 
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vector field that satisfies = —dKg. We have 

u{Ut),fs{t))=u{Xo,fs{t)) (3.106) 

Aiiit), m) = A{Xo, fit)), (3.107) 


which implies that the time derivative of the vector potential simplifies to 

^^Aiiit), fit)) = -AiXo, %)/(t)). (3.108) 

Thus, the time-dependent Hamilton equation simplifies to 

i^4i)d + -^7r*AiX,, fit))) - -^7r*AiX,, %)/(t)) + d— = 0, (3.109) 

which is merely a relationship between ^(t) and fit). 

The relationship between ^(t) and fit) can be resolved explicitly as follows. By applying fiber 
integrals to the Vlasov equation, /^(t) = —L^^(^t)fsit), we obtain the fluid equation for the particle 
flux Usit) = ui^sit), fsit)), 

^ * Usit) = -diviTsit))^ - —i*ng)-Ait) - -^(*n^(t))d0(f) + — * i*Ugit) A *dA(f)), 
dt rrig c rris nig 

(3.110) 

where 


risit) = nifgit)) (3.111) 

Ait) = Aiiit) Jit)) (3.112) 

0(t) = $(/>)) (3.113) 

Tsit) = Tifgit)) 


(3.114) 
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and T is the stress tensor operator. Given a pair of 1-forms a, (3 & the defining relation for 

the stress tensor operator is 


T(/)(«,/3) = * 




(3.115) 


By multiplying the fluid equation for Ugit) by and then summing over species, we then obtain 


4* j«(t),/(i)) = 

Ns . .. Ns 2 

- ^e,div(T,(f))‘' - A *dA(f)) 

S=1 S = 1 ^ ^ 

(3.116) 


as an equation for the time-derivative of the current density. Here we have introduced the local 
plasma frequency = Up{f{t)), where 


Ns , 2 

S=1 


m.o 




(3.117) 


Finally, by applying the operator —^GUt to both sides of Eq. (3.116), we obtain a linear operator 
Ej : V* ^ fl^(Q) that gives the inductive electric field in terms of the distribution function. 


ElU) 


1 -1 


1 + ^GUrUjlif) 


^GUTul{f)[Oif)], 


(3.118) 


where O : V* ^ (Q) is given by 

0{f) = d4(/) + V (e.diY(T(/))‘ - A , /.) A /)]) . (3.119) 
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The operator El is preeisely what is neeessary to express the Eulerian phase spaee veloeity 
in terms of the distribution funetion. Indeed, we have 

Ut) = ^s(/(t))-' - es7r*Eiifit))^ , (3.120) 

where, for eaeh / e V*, Us{f) is the sympleetie form on TQ given by 

= -d (ras'd + ^7r*A{Xo, /)j . (3.121) 

Note that a;s(/) = —d9s{f), where 

9s{f) = ras'd + ^7r*^(X„, /). (3.122) 

With the relationship between $,s{t) and f{t) resolved, we ean now substitute it into the AVD 
Euler-Lagrange equations. When this substitution is performed, the AVD equations beeome a 
first-order ODE in the variables (g(t), x(t), fo(t)) e G x V x V*. Explieitly, we have 

- es7r*Ei^ o g^(t) (3.123) 

Xs(t) = gs(tr(^9s(X,)-^^ (3.124) 

fo,s(t) = 0, (3.125) 

where the operators 9s, Ej, and the funetional derivatives SH/dfs are evaluated at g(f)*/o(f) G 
V*. In writing the equation for Xs(t), we have made use of the faet that ~ is 

a seeond order veetor field. 

As is true of first-order ODEs in general, this first-order ODE for the variables (g(t), A(t), fo{t)) ^ 
G X E X E* is identifiable with a veetor field E on G x E x V*. Setting Z{t) = (g(t), x{t), /o(f)), 


gs(f) = tu. 


-1 

5 


6H 

dfs 
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V is defined by the relation 


^Z(t) = V(Z(t)). (3.126) 

dt 

Because the AVD equations can be written in this form, it follows that the submanifold Po C 
T{G xV X V*) that serves as the AVD phase space is given by 

Po = {{Z, Z) e T{G xV xV*)\Z = Y{Z)} ^GxVxV*. (3.127) 

3.3.4 Step 4: Derivation of boundary symplectic form on Vo 

We will denote points in G x V x V* with the letter Z, i.e. Z = {g,x, /o). The vector field Y 
defines a time-independent flow map ■ Po = G x V x V* ^ Po = G x V x V*, which is 
characterized by the relations 


)?o = idp„ (3.128) 

^MZ) = Y{Z). (3.129) 

dr 

We can therefore define a mapping Sol : Po %i{G x V x V*), given by 

Sol{Z){t)=dt-t,{Z), (3.130) 

that sends initial conditions (at t = tp to their corresponding solution path in ^(G x V x V*). 

The mapping Sol can be used to pull back the augmented action functional S to the augmented 
phase space Po, thereby defining the restricted augmented action 


Sr, = Sor5. 


(3.131) 
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By examining the exterior derivative of the restrieted augmented aetion, we ean identify a sym- 
pleetie form, and therefore Poisson brackets, for the AVD equations. 

The exterior derivative of S-p^ is proportional to the free-endpoint variation of the augmented 
action functional S. Because this variation will be evaluated at a curve mG xV xV* that satisfies 
the AVD Euler-Lagrange equations, only the endpoint contributions to the free-endpoint variation 
will appear. Specifically we have 


dSp^ = ru_.e-e, 


(3.132) 


where 0 is the 1-form on Vo given by 

Ns . 

Q{Z)[5Z] = Y, {OsCgJo) ® {9s)Jo,s) ■ Vs + Xs5fo,s. (3.133) 

S=1 -^^<3 

Here, ps = Sgs o If we now differentiate Eq. (3.132) in the variable t 2 , keeping in mind that 
depends on t 2 via the upper limit of time integration, we obtain 

^Ly0, (3.134) 

where is the augmented Eagrangian pulled back to the augmented phase space via Y : Vo ^ 
T{G xV X V*),i.e. 

CpJZ) = {Y-C)(Z) 

= CiY(Z)) 

= X]y Ks{gs)*fo,s + — (^A{Xo,gJo),J{Xo,gJo)'j - -(^p{gJo),^{gJo)Y 


(3.135) 
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By applying the identity Ly0 = iyd© + diy0, Eq. (3.134) can be re-written as 


iyd© = —d'H'p^, 


(3.136) 


where 


= {Q{Y) - Cr.){Z) 



A{Xo,gJo),J{Xo,gJo)) + 7 ^ {p{g*fo),‘^{g*fo) 


1 

2 


(3.137) 


is the augmented Hamiltonian functional. Because the 2-form —d© turns out to be non-degenerate, 
Eq. (3.136) shows that the AVD equations are an infinite-dimensional Hamiltonian system with 
symplectic form — d© and Hamiltonian 

We will now calculate the Poisson bracket defined by the symplectic form —d0. We will 
proceed in two steps. Eirst, we will find an expression for an arbitrary Hamiltonian vector field on 
Vo by solving the equation 


iygd© = -dG, 


(3.138) 


for Yq given an arbitrary functional G : Po —^ Next we will identify an explicit expression for 

the Poisson bracket associated with —d© using the formula 


[F.GK = Lk,F, 


(3.139) 


Here, [•, -j-p^ is the augmented system’s Poisson bracket. Because calculating the Eie derivative in 
the last expression is simple, all of the nontrivial work will be done in the first step. 



Chapter 3. Boundary terms and Poisson brackets 


50 


A helpful tool for solving Eq. (3.138) is the phase spaee variational prineiple. This variational 
principle states that (fixed-endpoint) variations of the action functional &q : ^(Po) given by 

^ r ( 

are zero if and only if Z is a solution of the equation 

^Z{t) = Yc{Z{t)). (3.141) 

dt 

Thus, we know a priori that the Euler-Lagrange equations associated with the action ©g are 
Eq. (3.138). On the other hand, we can compute the fixed-endpoint variation of ©g directly, giving 


0[Z(f)] - G(Z(f)) dt 


(3.140) 


5©g(Z)[5Z] = 

|•t2 r 
Ns i>t2 r 

5=1 Jtq 

Ns^ /.t2 

5_i J t-i Jtq 


d{-7r*Am.fmxo)]-^^es- 


fsit) 


6G' 

ds i^sCqd^s T 
5gs. 


fs{t) ■ 775(f) dt 


Ssit)* Osi^sit)) + —Tr*A{^{t), fit)){Xo) - Xs{t) - 


SG 

Sfo,s 


Sfo,s{t) dt 


ut) - 


dXs{t) dt, 


(3.142) 


where / 5 (f) = gsit)*fo,si^), Vsit) = 5g{t) o g(f) \ and 0, is evaluated at /(f). By the phase 
space variational principle, if we set this variation equal to zero and then solve for Z(f) in terms 
of the functional derivatives of G, the result will be the solution to Eq. (3.138). We now turn to 
performing this task. 
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Just as when we formulated the AVD equations as a first-order ODE, solving for Z{t) involves 
dealing with an implicit (linear) equation for Indeed, the Euler-Lagrange equation given by 

0 = d /(i))(^o)) - ^9. - - C«)d9, - (3,143) 

is an implicit equation for because ^(f) appears in the first, second, and fourth terms. We 
will solve this implicit equation by working in the tangent lift of an arbitrary coordinate system g* 
on Q. Eirst we will express /(f)), which appears in the first term in Eq. (3.143), in terms 

of Z{t). Then we will express A(Xo, L^(j^/(f)), which appears in the second term of Eq. (3.143), 
in terms of Z{t). 

Solving for f{t)) — Eet gij denote the components of the metric tensor in our cho¬ 

sen coordinate system. The determinant of this matrix of components will be denoted \g\. 
Set 


Note that 


fs{t) = Fs-d{} A d?) A d-i) 
6 


Aiiit),f{t)) = Aid/ 


— = Qi dg* -f Qi dq\ 
ogs 


(3.144) 

(3.145) 

(3.146) 

(3.147) 


-d'i) A d-i) A d-i) = |g| dg^ A dg^ A dg^ A d/ A d/ A d/. (3.148) 

6 

In order to compute A.(^(f), /(f)), we will draw upon coordinate expressions for the current 
density operator J(^(f), /(f)) = /s(f))- As is readily verified, the fiber integrals 
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that give the partiele flux 1-forms, *u, ean be written as 

dg^, (3.149) 

where dq = dq^dq^dq^ denotes the “bare” measure on q spaee. The eurrent density 1-form is 
therefore given by 

Ns 

S = 1 

We will also draw upon the expression for ui that is implied by the dg*-component of Eq. (3.143), 
namely 


Fsulgij^/\i\dqj dq^ 


(3.150) 




Fsulgij^/\i\dqj 


gijul = -—-Ai + —Qi. (3.151) 

nis c nis 

By definition, the 1-form f{t)) is given by 

^(l(^),/W) = (/ FsKgij^/\i\dq^ dg^^, (3.152) 

where we have used the fiber integral identity given above. Both sides of this equation depend on 
^{t). However, upon inserting the expression for m* given above, we obtain 


1 / \ 4 /Ns 

Mm, /(*)) = --GnTC,M(«(«),/(i)) + —GYIt[ Ui)) 


J C 


S =1 


(3.153) 


= -^ght MMm,f(t)) + Amm/sg], f(t)), 


(3.154) 
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where u \5Q/5g] is the element of 0 given by (ce and we have used 

the identity 





(3.155) 


We therefore arrive at the simple conclusion 


-1 


Aim Jit)) = 1 + fit)) 

= n^iht)) \Am[SG/5g], fit))] , 


(3.156) 

(3.157) 


where we have introduced the squared refractive index operator 


n\f) 


1 + -^GIlT(2!pif) 


(3.158) 


Solving for AJo, L^^^-^ff)) — In order to solve for AJo, L^^^-^ff)) in terms of Zf), we 
will make use of coordinate expressions for the 2 -form d 6*5 and the 1 -form *uiXo, L^^(^t)fsit))- 
We have 


de. 


*uJo,Lut)fsit)) 


1 6 

rusgij df A dg-^ -f ( rUsTjkJ -f ) dg'" A dg^ 



(3.159) 

(3.160) 


where 


(3.161) 
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are the Christoffel symbols of the first kind, 

^Bkjdq'^ A dg-’ = dA{Xo, f{t)), 
and 

r“= j uU'F,^\dq. 

We will also employ a coordinate expression for the inverse of the 2-form Us = —dOg. 
1-form on TQ, a = aj dq^ + dj dq\ 


where 



} dq'- 


rri^c 


a 


+ 


rug^ ) dq^ ’ 


cm 


Blj — - q^i^jki — Tfeji] -f Bkj. 


Set A{t) = A{Xo, By definition, the 1-form A{t) is given by 


471 


Ns 


A(t) = —GUt Y.^g*u{Xo,L^A)Ut)) 


. 5 = 1 


If we decompose $,g{t) as 


- {A{t )), 

c 


(3.162) 

(3.163) 

Given a 

(3.164) 

(3.165) 

(3.166) 


isit) = Vg{t) 


(3.167) 
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where Usit) is the veetor field on phase spaee defined by 




6G gs{t)^{6G/6xs) a , * 


. ... Ss + —t^*A{Xo, g(t)*(5G/5x)) 
fsit) c 


- d (-^7r*n2(/(t)) A{u: ^[5G/5g\, f{t)) (X^)) , 


(3.168) 


then we ean deeompose A{Xo, into two pieees, 

A(t) = ^GUt e, * u(X„ j - GIIt 1 ^ * n{X„, 1 . 

(3.169) 

The first term on the right-hand-side of this expression is given entirely in terms of Z (t) beeause 
the veetor field Ug (t) only depends on Z (t). The seeond term on the right-hand-side involves the 
quantity we are trying to solve for, A{t). Using the eoordinate identities given earlier, the 1-form 
*u{Xo, can be expressed as 


777-c 


(3.170) 


Therefore Eq. (3.169) simplifies to 


Ns 


47r / \ 1 

A{t) = -^GUt \^^eg*u{Xo,L^xt)fs{t))j - -^GUt (a;|^(f)) , 


(3.171) 


whieh provides us with the expression for A{t) = A{Xo, we have sought after: 


AiX,,L^^Jit)) = n\fit)) AiXo,L^it)fit)) 


(3.172) 









Chapter 3. Boundary terms and Poisson brackets 


56 


With these expressions for and A{^{t),f{t)) in hand, we can finally write 

down an explicit expression for Yq. Set Y{g, x, fo) = {g, X, fo, 9 g, Xg, (/o)g)- We have 


{9s)g°9s ^ ^ (^*n^{9*fo) A{Xo, 


(3.173) 


(Xs)g = 


SG 


Sfo,s 


+ 9*s [Os { i^s { +^'n-*n‘^{9*fo) A{u ^[6G/6g], gjo) (Xo) 


if )g = ^ 

\JO,SjG r- ? 


(3.174) 

(3.175) 


where the operator z/^ : C* x g* x C* —)■ X(TQ) is given by 


Us{f;a;h) {as + Y0s + ^TT*AiXo,h)-d(^^7i*n^{f{t)) Aiuj ^[a\,fit)) (Xo) 


(3.176) 


It follows from these expressions together with Eq. (3.139) that the Poisson bracket on the 
augmented Vlasov-Darwin phase space can be written down immediately. However, the most 
obvious form of the bracket is not manifestly antisymmetric. After some toil, an antisymmetric 
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expression for the bracket can be found. It is given by 


[f,gk(z) = 



where iis : 0 x C* —)■ C°°(TQ) is a non-linear operator given by 


Uif) 




andcu 


G 0 is given by 


u 


-1 




(3.178) 


(3.179) 


In deriving this expression for the augmented system’s Poisson bracket, we have made use of the 
fact that the operator n^{f)GIlT is self-adjoint. See the next brief subsection for a proof of the 
self-adjoint property. 
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Properties of the squared refractive index operator 

The squared refractive index operator is defined by 


^^(/) = 


n -1 


1 + -GUtuIU) 


( 3 . 180 ) 


where a>p(/) is the operator that simply multiplies by the (squared) local plasma frequency. n^(/) 
operates on differential forms over Q. This definition may be perplexing because it is not immedi¬ 
ately obvious that A = l + jzGIlT(^p{f) should be an invertible operator. To see that it is, suppose 
that A has a non-trivial null eigenvector ao- Then ao would have to satisfy 


^GUTuliDao = -ao. ( 3 . 181 ) 

In particular, Uo would have to be an eigenvector of the operator ^Gnra)p(/) with eigenvalue —1. 
This is impossible for the following reason. Define the weighted inner product 


{a,P)^= / {aA*/3)ul{f). (3.182) 

JQ 

Because the operator GAt is the product of non-negative definite operators, it is non-negative 
definite itself. Therefore, for each a, 

^GAtuICDo^ = GATulCf)a) > 0. (3.183) 

But this contradicts our assumption that there is a non-trivial null eigenvector of ao- Indeed, 

ao, ^Gnra)p(/)ao\ = - (ho, To) < 0. (3.184) 

^ /a; 

It follows that the operator A is invertible, and that n^(/) is well-defined. 




Chapter 3. Boundary terms and Poisson brackets 


59 


The most basic property of is that it commutes with E = To see this, set 

C = n‘^if)E - En^if). We have 

(1 + E)C = E-{1 + E)En\f) 

^(1 + E)C{1 + E) = E{1 + E)-{1 + E)E = 0 

^C = 0. (3.185) 

While n^{f) is not self-adjoint as on operator on the standard space of differential forms, 
it is self-adjoint as an operator on the weighted space defined by the weighted inner product 
introduced earlier. To see this, let a and (3 be arbitrary /c-forms in Set ai = 'n3{f)a and 
/Si = n^(/)/9. We have 

(a,n=(/)/?)_^= ^ l + iGnTW=(/) ai,/3,^ 

= (£,^(/>i,/?i) + {^iGnTwJ(/)ai>J(/)/3i\ 

= (£,2(/)ai,/?i) + ^wj(/)a„ iGnTwJ(/)/3i\ 

= /i.J(/)ai, l + iGnTW=(/) ft) 

—/ai, 1 +—GUruj^if) 

\ L J / oj 

= (^n\f)a,/3)^, (3.186) 


where we have used the fact that GUt is a self-adjoint operator on E^. 




Chapter 3. Boundary terms and Poisson brackets 


60 


When we combine the last two properties of the squared refractive index operator, we obtain 
the important result that n‘^{f)GIlT is self-adjoint on Indeed, 


a,n\f)GUTf3) = 


a 


^pif) 

= { n\f) 


,n\f)GTlTP 


a — ^ 


-^GliTU^Jf)- 




-p\ 
. 2 / 


a 


/S 


= Igut^IUWU)- 

\ Wp(/) ^pif) 


= (n\f)GUTa,/3 


(3.187) 


3.3.5 Step 5: perform Poisson reduction to obtain bracket on physical phase 
space 

We have now identified the Hamiltonian, 


'HrSa^xJo) = I i9s)Jo,s + — (^A{Xo,gJo),J{Xo,gJo)'j + - (^p{g*fo),^{gJo)'j , 

(3.188) 


S = 1 


and Poisson bracket for the augmented Vlasov-Darwin system. The Jacobi identity is satisfied 
because the bracket has been obtained by inverting the symplectic form —dS. The modifier “aug¬ 
mented” is appropriate because the dynamical variable x has no direct physical meaning. On the 
other hand, the variables g and /« together comprise an element of the physical Lagrangian (as 
opposed to Eulerian) phase space; g gives the configuration of particles in the single-particle phase 
space TQ and /« gives the reference phase space density. The purpose of the additional variable x 
is to extend the Lagrangian phase space just enough to allow for a non-degenerate Poisson bracket. 
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The appearance of the variable x perhaps seems awkward at this stage. However, observe the 
following. The set C is a Lie group under addition that is a symmetry group for the AVD equations. 
Specifically, for each 5x ^ L, we can define a mapping ■ Vo ^ Vo given by 

Tsxig, X, fo) = {9,X + Sx, fo), (3.189) 

which clearly satisfies the defining properties of a group action, 

^5x1+<5x2 = ^<5X1 ° ^<5X2 (3.190) 

To = idp„. (3.191) 

This group action leaves the augmented Hamiltonian and the augmented Poisson bracket 
[•, ■]'p^ invariant in the sense that 


T*s^nv^=nv^ (3.192) 

TjilF.Gip. = |r;jF,T;,G|p., (3.193) 

for arbitrary 5x ^ V and functionals F, G on the augmented phase space. These properties are 
quick to verify. Therefore we can define a Hamiltonian and Poisson bracket on the first reduced 
phase space Vl = T’o/L, i.e we can perform Poisson reduction. 

The details of applying Poisson reduction to pass from the AVD phase space to the Lagrangian 
phase space Vl follow. The Lagrangian Hamiltonian, —)■ M, is uniquely determined by 

requiring 


X, fo)) = Vp^ig, X, fo) 


(3.194) 
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for each {g, x, fo) ^ Vo- Here the projection map til ■ Vo ^ Vl is given by 

T^L{g,xJo) = {g,fo)- (3.195) 

The Lagrangian Poisson bracket, [•, is defined by requiring 

7rl[d,&]r, = Kd,7rliSW (3.196) 

These definitions make sense because the augmented Hamiltonian and Poisson bracket are invari¬ 
ant under the the action of V. We find that the Lagrangian Hamiltonian is given by 

Vv^igJo) = / ^sgsjos + — (^A{Xo,gJo),J{Xo,gJo)'j + - (^p{gJo),^{gJo)'j ■ 

s=l P 

(3.197) 


The Lagrangian Poisson bracket is given by 

Ns . 

l5.ek = E/ 

s=l JTQ 

(3.198) 


Sgs 




I > 9* Jo , ^ 

og J J ogs 


dii. 


-1 AA ~ 

^ [^] ^9* fc 

\^9 J 


gs*fo. 


This bracket and Hamiltonian give a Hamiltonian formulation of the Vlasov-Darwin equations in 
Lagrangian labeling. 

Now we will pass from Lagrangian labeling to Eulerian labeling by applying Poisson reduction 
a second time. The set G is a symmetry group of the Vlasov-Darwin system in Lagrangian labeling. 
Specifically, for each h e G, we can define a mapping Vl given by 


VhivJo) = {goh, h*fo), 


(3.199) 
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that satisfies the defining properties of a right group aetion, namely 

^hioh2 ^ ^h2 ° ^hi (3.200) 

f?idp=idp^. (3.201) 

The Lagrangian Hamiltonian and Poisson braeket are eaeh invariant under this group aetion, in the 
sense that 


RlHr, = Hr, (3.202) 

Rl[J^:0]r, = [RlJ^,RlGW (3.203) 

The invarianee of the Lagrangian Hamiltonian is quiek to verify. The invarianee of the Poisson 
bracket follows from the identity 

(9. /•) = ° (S'ZO-*) 

Therefore, the Lagrangian Hamiltonian and Poisson bracket define corresponding quantities on the 
Eulerian phase space, Ve = Vl/G ’^V*. The Eulerian Hamiltonian is given by 

Ns /* I 3 

= / K,f. + -{A(X„f),J(X„f))+-(p(f),i(f)). (3.205) 

s=l 

The Eulerian Poisson bracket is defined in terms of the Eulerian projection map, tte ■ Vl ^ Re, 
given by 


^e{ 9, f o) 9*fo, 


(3.206) 
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which satisfies the important property tte o Rh = We have for functionals ^ —)■ M, 

= E ^ (-V/./. /) . f - /) }/.. (3.207) 

where g and /o are any group element and 6-form that satisfy / = g^fo, {•, - js is the Poisson 
bracket associated with the symplectic form Ug, and G } is given by 

(AV/tf), = = ‘^;‘(d3A3/.). (3.208) 

3.4 The gyrokinetic Vlasov-Maxwell system 

The gyrokinetic Maxwell-Vlasov system is most naturally defined by specifying its action, which 
is the sum of the net gyrocenter action and the Maxwell action. If a gauge-invariant form of the 
gyrocenter transformation is employed, we have 

SjSz.a,4>) = 

(3.209) 

where ^{t) = g o is the multi-species Eulerian phase space fluid velocity, f{t) = g(f)*/o is the 
multi-species gyrocenter phase space density (a collection of 6-forms), E{t) = —d 4 >{t) — A{t )/c 
is the electric field 1-form, B{t) = dA{t) is the magnetic field 2-form and the various Lagrangian 
functions are defined as follows. Cp is the “free gyrocenter” Lagrangian (in analogy with the notion 
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of a free-particle Lagrangian) given by 


S=1 





(3.210) 


where is the guiding center 1-form and /C* is the gyrocenter kinetic energy functional 


ICs{f,E,B)= [ fK,{E,B). (3.211) 

Jtq 

Here Kg : fl^(Q) x df2^((5) —)■ C°°{TQ) is the gyrocenter kinetic energy function, which has 
a non-local dependence on the electric and magnetic fields. £int is the “interaction” Lagrangian 
given by 


Ns 

s=l 


/ fs 

—'k*A{Q - egTT*(p 

Jtq 

1 c J 


(3.212) 


where the current density 1-form and charge density operators, J and p, were introduced in the 
previous section. Finally, Cuax is the free electromagnetic field action 


Cua.{E,B) = — {E,E)- — {B,B). 
8n 8n 


(3.213) 


Explicit expressions for H and Kg will not be necessary; it is enough to know they can be found in 
principle. 
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The Euler-Lagrange equations associated with the gyrokinetic action are given by 


d \ f ^ -^int) 

* + 

1 d 6{Cp + CMax) + -^Max) 

cdt Je Jb 

5(^C,p T Cjviax) 


= d 


6{Cp + Cint) 


5fs 


fs 




5A 

SCint 


5E 


Scj) 


Upon calculating the relevant functional derivatives, these equations reduce to 


i^^qdSf = -Sf - dUf 


~D(t) - SH(t) = -Ej(|(t),/(i)) 

c dr c 

5D{t) = -47rp(/(t)), 


where we have introduced the gyrocenter 1-form and Hamiltonian, 


Ef = H, + 


Uf = Ks + e^7r*0, 


the 1-form 


D(t) = E(t) - iJE = £,(/((), E(t), B(,t)), 


and the 2-form 


H{t) = B(t) + iJE = B{t)). 


(3.214) 

(3.215) 

(3.216) 

(3.217) 

(3.218) 

(3.219) 

(3.220) 

(3.221) 

(3.222) 


(3.223) 
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Note that in these expressions E(t) = —d<p{t) — A(t)/c, B{t) = dA{t), and 1C = 
which implies dE(t) = —B{t)/c. Equation (3.217) gives the gyrocenter equations of motion in 
Hamiltonian form. Equation (3.218) is the “macroscopic” Ampere equation (written in terms of 
forms). Einally, Eq. (3.219) is the macroscopic Gauss equation. 

The reference gyrocenter phase space density fo is not a dynamical variable in this formulation 
of gyrokinetics. However, we can elevate fo to the status of a dynamical variable by embedding 
gyrokinetics within a slightly larger system. In particular, if we define the augmented gyrokinetic 
action 


(f),fo,x) 


N, 




S=1 



fos{t)Xs{t), 


the resulting Euler-Eagrange equations are given by 


(3.224) 


= -Sf - dWf 

(3.225) 

c 

(3.226) 

6D{t) = -47rp(/(f)), 

(3.227) 

fosif) = 0 

(3.228) 


(3.229) 


The first four equations in this set decouple from the fifth and reproduce the gyrokinetic Vlasov- 
Maxwell dynamics. The fifth equation defines the dynamics of the additional variable x(f). We 
will refer to this larger system of equations as the augmented gyrokinetic Vlasov-Maxwell equa¬ 
tions, or the AGVM equations for short. 

It may seem awkward to introduce the additional dynamical variable x(f). However, x{t) 
will be a help rather than a hinderance as we continue our derivation of the gyrokinetic bracket. 
Moreover, the evolution equation for x{t) is very much reminiscent of the evolution equation 
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for the phase of the quantum wave funetion in the WKB approximation. I don’t think this is a 
coincidence; a semiclassical treatment of gyrokinetics would most likely produce an equation very 
similar to Eq. (3.229). This is strongly suggested by considering the gauge invariance properties of 
the augmented gyrokinetic action. Because fos{i) depends on time, under the gauge transformation 
A{t) A{t) + —)■ — <p{t)/c the augmented gyrokinetic action is not weakly 

invariant (it changes by more than temporal boundary terms). However, if we also change Xs{t) 
according to Xs{t) —^ Xs{l) — the gyrokinetic action is left weakly invariant. 

Thus, Xs{t) changes in the same way as the quantum phase under a gauge transformation. 

The AGMV equations determine an evolution equation for the variable Z{t) = (g(f), A{t)^ 
D{t), /o(f), xil)) that is first order in time. To see this, first choose the gauge 4>{t) = 0. In this 
gauge, B{t) = dA(f) and E{t) = —A{t)/c. Therefore, if we introduce the electric field operator 
E, which is defined by the implicit equation 

D = b{f,E{f,D,B),B), (3.230) 


we can write 


AA(f) = -cEifit), D{t), dA{t)), (3.231) 

where f{t) = g(f)*/o(f), which gives the time derivative of A(t) in terms of Z{t). For the time 
derivative of g(f), we note that 


Ss{t) = (3.232) 

and that is a functional of E and B according to Hamilton’s equations, i.e. ^s(f) = 
Xf^{E{f{t),D{t),dA{t)),dA{t)) for a functional : H^(Q) x H^(Q) —)■ X{TQ). Because 
^s{t) can be expressed in terms of Z{t), it follows that the macroscopic Ampere equation gives 
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the time derivative of D{t) in terms of Z{t). For the same reason, the time derivative of Xs{i) can 
be written in terms of Z{t). Finally, the time derivative of /os(t) is trivially a funetional of Z{t). 
It follows that the variable Z{t) obeys an (infinite-dimensional) autonomous first order ODE 

Z{t) = Y{Z{t)), (3.233) 

where y is a veetor field onVo = G x f2^((5) x f2^((5) x V* x V. We will refer to the spaee Vo 
as the AGVM phase spaee. We will use ■ Vo ^ Vo to denote the (formal) flow map assoeiated 
with Y. 

By the existenee and uniqueness of solutions to first order ODEs, for eaeh point Z G Vo, there 
is a unique path Z sueh that Z{ti) = Z and Z{t) = Y{Z{t)). In terms of the flow map this 
path is given by 

Z{t)=:St-tAZ)=^o\{Z){t), (3.234) 

where we have defined the funetion Sol : Vo —)■ ‘^{Vo)- The augmented gyrokinetie system 
therefore gives us a natural way of mapping points in the augmented phase spaee into a path 
spaee. But reeall that the augmented gyrokinetie aetion maps points in a path spaee into the real 
numbers. This suggests that we ean eonstruet a speeial real-valued funetion on Vo by eomposing 
the gyrokinetie aetion with the funetion Sol in some sense. Aetually, eomposition doesn’t make 
literal sense because the augmented gyrokinetie action is a functional defined on a path space 
that differs from ‘^{Vo)', the argument of the augmented gyrokinetie action is a path of the form 
(g, A, (p, fo, x) ^ ^ X V* X V). Nevertheless, there is a simple mapping 

n ; ^{Vo) 'P(G X Di(g) X D0(g) xV*xV) given by 


n(g, A, D, fo, x) = (g, A, 0, fo, x), 


(3.235) 
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which is motivated by our choice of the gauge 0(t) = 0. We can still therefore define a special 
real-valued function on Vo 


^P„(Z) = ^(n(Sol(Z))). (3.236) 

We will refer to S-p^ as the restricted augmented action. Note that implicitly depends on ti and 
^2 as these appear in the action integral’s limits of integration. 

It is extremely interesting to study the (infinite-dimensional) exterior derivative of the restricted 
augmented action. A simple direct calculation shows 

dSvXSZ) = - dr,, (3.237) 

where 0-p^ is a 1-form on Vo given by 

QvX^Z) = ^ 0 g,Jos) • Vs + XsSfos^ - ^ {D, 6A) . (3.238) 

Here v = 5gog~^. This identity becomes very interesting indeed when both sides are differentiated 
with respect to ^ 2 - We have 

= 5;._,.(iyder,. + d(e„,(K))), (3.239) 

where Sp^ is given by 

Sp^{Z) = gjo, E, dA) + 0) + Cu..{E, dA). (3.240) 

Here we are using the short hand notation = X^^{E{g^,fo, D,dA),dA) and E = 

E{E{g^fo, D,dA). Equation (3.239) tells us that the augmented gyrokinetic dynamical vec- 
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tor field V obeys an infinite-dimensional version of Hamilton’s equations, 


(3.241) 


where the elosed (and non-degenerate, as it turns out) 2-form = -dO-p^ and the energy fune- 
tional = 0p^(V) — Sp^. It follows that the AGVM equations have a Poisson formulation on 
Vo with a Poisson braeket given by inverting the sympleetie form H and a Hamiltonian funetional 
given hy 8p^. 

In order to derive an explicit expression for the AGVM Poisson bracket [•> -]vo , the simplest 
approach is to first derive an expression for a general Hamiltonian vector field, Xq, using the 
infinite-dimensional phase space variational principle 

and then calculate the Poisson bracket using the formula 


Z{t)) - G(Z(t)) Mt = 0, 


(3.242) 


[F,GK = Lx,F. 


We find that Xg = (fc o 9, ^g, Dq, (/o)g, Xg), where 




.4= = 

(5G 

= 47rc— - 47rJ(fG,^*/o) 


5A 


{fos)G — 


6G 

SXs 


ixs)G = 9 mms)G])- 


6G 

Sfos 


(3.243) 


(3.244) 

(3.245) 

(3.246) 

(3.247) 


(3.248) 
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Therefore the AGVM Poisson braeket is given by 


5F ^ 6G ^ JG . ^ 

- -47re,7r —, -47re,7r — g^Jo 

dgs oD dgs oD 


1f.cK = E/ 

s=l JTQ 


N, 


E 


(5F (5G 5G (5F 


'■^JtqS/osSxs SfosSxs) Ww^SA/ \SD^SA 


^ ’ 6gs * 6D 


(3.249) 


where Bf^ is the gyroeenter Poisson tensor. 

The AGVM Poisson braeket in eonjunetion with the energy funetional 


Sr.iZ) = (er^y) - SrJ{Z) 

= fCCgJo, E, dA) + (P, e) + ^ (p, e) + ^ {dA, dA), (3.250) 

east the AGVM equations as an infinite dimensional Poisson dynamieal system. Here we have 
introdueed the polarization operator 




1 

An 


{D-E{f,D,B)). 


(3.251) 


Note that the AGVM Poisson braeket is derived by inverting a sympleetie form, and so it does not 
have easimirs. The main drawbaek of this Poisson formulation of gyrokineties is that it eontains 
a superfluous dynamieal variable x. A seeondary drawbaek is that it uses the gauge-dependent 
veetor potential A as a dynamieal variable. A tertiary drawbaek is that it is expressed in Lagrangian 
labeling. 

All of the drawbaeks of the augmented Poisson formulation of gyrokineties ean be removed 
using Poisson reduetion. To refresh the reader’s memory, the basie premise of Poisson reduetion 
is that if there is a Lie group H that aets freely on a Poisson manifold P in a manner that leaves 
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the Poisson bracket and Hamiltonian function invariant, then the bracket and Hamiltonian induce 
corresponding objects on the quotient P/H. The bracket on P/H is known as the reduced bracket 
and the Hamiltonian on P/H is known as the reduced Hamiltonian. In our case, the poisson 
manifold will he P = Po equipped with the AGVM bracket, and the Hamiltonian function will be 
The Lie group will be 77 = 72° (Q) x {G xi V). The first factor in H corresponds to gauge 
symmetry, the second to particle-relabeling symmetry, and the third to symmetry with respect 
to translations in the augmented variable x- The quotient P/H will be the gyrokinetic Eulerian 
phase space, i.e. the space of triples (/, 72, B). The reduced bracket and reduced Hamiltonian on 
the Eulerian phase space will provide us with a Poisson formulation of the gyrokinetic Vlasov- 
Maxwell system without any of the drawbacks of the AGVM Poisson formulation. 

Instead of applying Poisson reduction to the entire symmetry group of the AGVM system 
77 = 72°(Q) X (G XI V), we will apply the Poisson reduction procedure three times, eliminating 
one factor of 77 in each step. The general theory behind breaking a symmetry group into subgroups 
and then applying step-wise Poisson reduction is known as Poisson reduction by stages. We will 
not need to draw upon any of the general results from this theory. Eirst we will factor out the 
translations in x, which are embodied by the additive Eie group V. This will give us a Poisson 
formulation of gyrokinetics on the space Vla = T’o/L, which consists of tuples {g, A, 72, fo). Pla 
is the gauge-dependent Eagrangian phase space for gyrokinetics. In the second reduction step, we 
will factor out the gauge symmetry, which is generated by a second additive Eie group 72° (Q). The 
result will be a Poisson formulation of gyrokinetics on the gauge-independent Eagrangian phase 
space Pl = Pla/^^{Q)- Einally, we will quotient by the particle relabeling symmetry group, 
G, which will provide us with a Poisson formulation of gyrokinetics on the Eulerian phase space 
Pe = Pl/G. 
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The Poisson bracket on the first reduced phase space, Vla = Vo/V, is given by 


1F,GK,(9,AA/o) = E / 

s=l JTQ 


[sg - SD'S^r W I 


, , ^ 6F 5G\ / 5G 5F 

*’’^^^sd'Ja)-{w’sa 


(3.252) 


The reduced Hamiltonian on Vla is given by 


{~g, A, D, Q = }C{~gJo, E,dA) + {P, E) + ^ (e, e) + ^ {dA, dA). (3.253) 


Note that this reduced bracket and Hamiltonian are both invariant under time-independent gauge 
transformations A —>■ A + dtp, where p E 

The Poisson bracket on the second reduced phase space, Vl = Vla/GE{Q), which is the 
gauge-independent Lagrangian phase space, is given by 


Ns . 

1f,gk(j,AB./o) = E / 

s=i Jtq 


f 6F ^ JF 5G ^ JG 


9s* fo 


Attc 


5D'^5B/ \5D'^5B 


(3.254) 


The Hamiltonian on is given by 


D, B, Q = ICCgJo, E, B) + {P, E) + ^ (^E, e) + ^ {B, B). 


(3.255) 


Note that the 2-form B is required to be exact (ignoring possible Homological complications). 
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Finally, the Poisson bracket on the Eulerian phase space, Ve = Vl/G is given by 


|F, G]„,(/, C, B) = ^ / Bf - 4xe, 

^ JtQ V ofs 


JG ^ JG\ , 


, f/6f JG\ /6G ^6f 
^ 5b) \5D'^5B 


(3.256) 


and the Eulerian Hamiltonian is given by 


SrSl D, B) = ICif, E, B) + (p, e) + ^ (e, e) + ^ {B, B). (3.257) 


Note that this bracket, and all of the previous brackets satisfy the Jacobi identity by construction. 


3.4.1 The augmented gyrokinetic Vlasov-Maxwell symmetry group 

The basic symmetry group of the AGVM equations is H = x (G x C), which as a set is 

simply 12°(Q) x G x C. Eor technical reasons, we will regard 12°(Q) as the functions on Q that 
vanish at infinity. Some basic properties of P are: 

• The group identity e = (0,idrQ,0). The group product of si = and si = 

( 992 ,^ 2 , 72 ) is given by 

<si * S 2 = {^1 + (p2, hi o /i 2 , h* 2 Ti + r 2 ). (3.258) 

The inverse of s = ( 99 , h, f) is given by 

s~^ = {—(f, h~^, —Kf). 


(3.259) 
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• i7’s Lie algebra f) = x g x V. We will denote typical elements of f) with the symbol 

X = (^(f, (, 5x)- The adjoint action of 77 on [) is given by 


Ad^a; = {5(p, KC, K[L^f + 5r]). 


(3.260) 


The Lie bracket is therefore 


[xi,X2] = (0,(3.261) 

• The dual to i7’s Lie algebra i)* = Q^{Q) x g* x V*. We will denote typical elements of ()* 
with the symbol /r = (p, a, /). The coadjoint action of H on [)* is given by 

Ad*/i = (p, /i*[q; — df ®/], 7*/). (3.262) 

The conditions for a functional A : [)* ^ M to be invariant under the coadjoint action are 
therefore 


L5A/Saf = 0 (3.263) 

6A 

LsA/s&a + d-^ 0 / = 0. (3.264) 

dj 

Note that these conditions can be regarded as first-order functional partial differential equa¬ 
tions for the Casimirs of the Lie Poisson bracket on ()*. 

• There is a right i7-action on Vo given by 

Rs{g, D, A, fo, x) = {g oh,D,A + dp, h*fo, h*{x + (e/c)p*7r» + r). (3.265) 
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The right infinitesimal generator is given by 

xr{Z) = {gX 0^,0, d5ip, L~Jo, L^x + {e/c)g*7^*S(p + 6f). (3.266) 

Note that the right infinitesimal generator satisfies B*gXR = (Adsa;)i{ for eaeh s G H. 

In order to identify the eonservation laws associated with this symmetry group, we will now 
consider the invariance properties of the 1-form A straightforward calculation shows 

+ dc„ (3.267) 


where 

Cs(Z^ \ '' I 71 gs*fos T hg^Tsfos- 

Differentiating this equation in s at e G 77, we obtain 


(3.268) 


ia;«d0 = -d (0(a;ij) - Sc{x)), 


(3.269) 


for each a; e f). Here we have introduced the mapping (ic : t) — )■ C°°(Vo) given by 

(3.210) 
(3.271) 

Q{xr) is given by 


Sc(x)(Z) = — Cexp(ex)(^) 


iVs p 

j ^Tgfos + - (^p{gJo), • 


TV, 

Q{xr){Z) = Y ^9*s^s^ ® fos - dx 0 fos) -Cs- — {SD, 6(p). (3.272) 

g^i Jtq 47 rc 




Chapter 3. Boundary terms and Poisson brackets 


78 


It follows that a momentum mapping m : Vo ^ i)* associated with the action Rs is given by 

D, A, fo, x)= (^ - ~ - dx] 0 fo, ■ (3.273) 

It is straightforward to verify that m satisfies the following equivariance property 

mo Rg = Ad*-i o m, (3.274) 

i.e. translating in Vo along the i7-action only changes the value of m by shifting it along a coadjoint 
orbit in [)*. Thus, m : Po —^ is a Poisson map when [)* is equipped with its Lie-Poisson bracket. 

We can now formally write down many of the Casimirs of the gyrokinetic Vlasov-Maxwell 
bracket in Eulerian labeling. Because m is a constant of motion for the AGVM system, any func¬ 
tion of m is also a constant of motion. In particular, if A : ^ M is an Ad*-invariant functional 

on g*, A o m is a constant of motion. The constant of motion Ao m satisfies /?* (A o m) = Aom 
by the equivariance of m and the Ad*-invariance of A. Therefore Aom descends to the quotient 
VojH = Ve, i-e- there is a functional Ca'-Ve^^ uniquely characterized by the formula 

= Aom, (3.275) 

where tte ■ Vo ^ VolH = Ve is the quotient map. Ca is a Casimir of the Eulerian gyrokinetic 
Vlasov-Maxwell bracket [•, because 


'^*e[^Ai'^]ve — [^*E^AiT^*e^]Vo 

(3.276) 

= [Aom,7r^F]p^ 

(3.277) 

= 0. 

(3.278) 
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3.5 The neutral Vlasov model 


The neutral Vlasov model was introdueed by Tronei and Camporeal in Tronei and Camporeale 
(2015). Its defining equations are given by 


dtfs + V ■ Vfs+—{E + vxB)- = 0 

rris 

dtB = -V X E 

V X B = fio'^es J vfs dv, J fsdv = 0. 


(3.279) 

(3.280) 

(3.281) 


We will first formulate the quasineutral model as an initial value problem without assuming 
V ■ = 0, whieh is not obviously implied by the Euler-Lagrange equations. We will work in 

terms of potentials in the temporal gauge = 0. I will also make use of a hodge deeomposition of 
the vector potential 


A = V X CK + VA, 


(3.282) 


where a. satisfies 


V ■ a = 0 (3.283) 

and A is defined modulo the addition of constants (i.e. A is really an equivalence class of functions 
whose elements differ from one another by constant functions). Given an A, there is a unique pair 
CK, A that satisfy the previous three equations. 

An equation for ct — According to the variational formulation of the neutral Vlasov model 
given in Tronei and Camporeale (2015), one of the Euler-Lagrange equations is p = 0, which 
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implies that V ■ J = 0. This is consistent with another one of the Euler-Lagrange equations, 

HoJ = V X B = —V X Aa. (3.284) 

Thus, we come to the conclusion that a. is uniquely determined by J (and therefore /) via 

a = /ioG[G[V X J]], (3.285) 

where G is the inverse of the Laplacian {G stands for Green operator). Note that A is not yet 
determined. 

An equation for dtX — If we take the time derivative of the Ampere equation, we obtain 

HoJ = -V X (V X E). (3.286) 

The time derivative of the current density can also be written in terms of the electromagnetic field 
using the momentum equation 

6 72 

dt{nsUs) + V ■Ts = ^{E + u,x B), (3.287) 

TUs 

where T is the stress tensor. 


We have 


T. 


vvfs dv. 


(3.288) 


j — —V • Q + SotOpE + eotOp { u ) x B, 


(3.289) 
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where 


Q = 


(3.290) 




{u) 



E 


n,et 


eoUis 


y a;2 

Z-^s ps 


(3.291) 

(3.292) 


If we now equate our two expressions for the time derivative of the eurrent density, we obtain an 
equation that funetionally relates the eleetrie field to the distribution funetion and the magnetie 
field, 


00 . 


1 

—V X (V X ^;) + ^; = -{u) X B + — 


)2 


(3.293) 


We will now use this expression to derive an expression for dtX. 
Note that the (/-dependent) linear operator E i—)■ BlE], where 


V[E] = —V X {V X E) + E, 


OJt. 


(3.294) 


does not have non-zero null eigenveetors. Indeed, if Eo is a null eigenveetor then we have 


0 < / '^\Eo\'^dx = 


j Eo ■ V X {V X Eo) dx = — J IV X Eol"^ dx < 0, (3.295) 


whieh implies that \Eo\ = 0 everywhere. Thus, we ean use the inverse of V to define an electric 
field functional £ = £{f, B), where 


£ 


V 


-1 


(n) X B 


eoOJt 


(3.296) 
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Because the vector potential determines the electric field according to the relation E = —dtA, the 
functional £ is related to dtX according to 

V-£{f,B) = -AdtX. (3.297) 

Because B can be regarded as a functional of / according to Eq. (3.285), we can now express dtX 
entirely in terms of /: 


dtX = -G[V-£{f,B)]. (3.298) 

The quasineutral model’s initial value formulation — We can now formulate the initial value 
problem for the quasineutral model. The state of the system is determined by the pair (/, A), where 
/ is the multi-species distribution function and A appears in the Hodge decomposition of the vector 
potential as above. The time evolution equation for (/, A) is given by 

dtfs = -V ■ V/, - —{£ + vxB)- V.fs (3.299) 

ms 

dtX = -G[V-£{f,B)], (3.300) 


where 


B = -fioG[V X J] (3.301) 

is expressed in terms of /. Note that the dynamics of / decouple from the dynamics of A. 

Remarks — This system of evolution equations admits solutions that do not satisfy the Ampere 
equation and the neutrality condition. However, using the definition of £, it is straightforward to 
verify that V ■ J is a constant of motion. Thus, if we choose initial conditions with V ■ J = 0, then 
this identity will hold for all subsequent times, and, as a consequence, the Ampere equation will be 



Chapter 3. Boundary terms and Poisson brackets 


83 


satisfied for all times (this second statement follows from applying the Laplacian to Eq. (3.301)). 
If, in addition to demanding our initial condition satisfies V ■ J = 0, we demand that our initial 
condition satisfies p = 0, then p will be zero for all times. Indeed, 

dtp = -V ■ J, (3.302) 

which implies that p is constant along solutions that initially satisfy V ■ J = 0. 

3.5.1 Phase space Lagrangian formulation of the neutral Vlasov model 

So far we have deduced that the phase space for the neutral Vlasov model is given by the space of 
pairs (/, A) subject to the holonomic constraint V ■ J = 0. Because we will need a slightly more 
precise notation in the discussion that follows, we will write the constraint as V ■ J(Xo, /) = 0, 
where the functional J(V, /) depends on the multi-species velocity field V and the multi-species 
distribution function / according to 

J{VJ) = Y,^s [ xjsdv, (3.303) 

S 

where Xg = Xs(x,v) is the spatial component of the phase space velocity Vg = (xgj'^s)- The 
vector field Xo is given by Xg = v and Us = 0, which implies 


as expected. 


J (^o, f) = eg j vfg dv, 


(3.304) 
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We will derive the Poisson braeket on this phase spaee in this seetion. We start by writing down 
a phase spaee Lagrangian for an augmented system, 

Ln{9,Fo,X,\9,Fo,xA) = Fo,X,'>^) -Fo, X, A) (3.305) 

Here g is a multi-speeies phase spaee fluid eonfiguration, i.e. an element of Diff(T(5); Fo is a 
multi-speeies referenee phase spaee density related to the referenee distribution funetion fo by 
Fos = fos where XI = dx A dy A dz A dvx A dvy A dvz', X is ^ multi-speeies funetion on 
phase spaee; and A is a funetion on Q modulo eonstant funetions. The multi-speeies phase spaee 
fluid veloeity is related to^byC = g o 9 ~^. The quantities g,Fo are eonstrained to satisfy 
V • J {Xo, F) = 0, where Fg = gs*Fos- The Hamiltonian funetional is given by 

n{g,Fo,x,X) = J2{Ks,Fg) + ^(V x At(Jo), V x At{Jo)) (3.306) 

where, Jo = J(3fo, F), At{Jo) is the transverse veetor potential given by 

At{Jo) = /ioG[G[V X (V X Jo)]], (3.307) 

and Kg = \mgV ■ v is the single-partiele kinetie energy. The Lagrange 1-form is given by 

0(L, Fo, X, A) = J2{MVs),Fg) + Fos). (3.308) 

5 S 

Here the singe-partiele Lagrange 1-form i9s is given by 


^g = {rrisV + es[A'r(Jo) -f VA]) ■ dx . 


(3.309) 
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The Hamiltonian vector field with Hamiltonian "H is determined by varying the action 
S-u = L-^ dt. More generally, the Hamiltonian vector field Xjr with Hamiltonian X is deter¬ 

mined by varying the action Sj^ = Ljr dt. The Poisson bracket [X, G] on the augmented phase 
space (i.e. {g, Fg, X: A)-space) is given by 


[F,G] = dF{Xg). (3.310) 

In order to derive the bracket, we will first find an expression for the general Hamiltonian vector 
field Xj: and then deduce the bracket from the previous formula. 

We will find Xjr by manipulating the Euler-Lagrange equations associated with the Lagrangian 
Ljr. In order to derive the Euler-Lagrange equations associated with the phase space Lagrangian 
Ljr, it is necessary to account for the constraint V ■ J(Xg, F) = 0. This holonomic constraint places 
a an awkward constraint on the variations of the phase space variables {g, Fg, x. A). I have found it 
very convenient to work with unconstrained variations and a Lagrange multiplier. That is, instead 
of varying the action Sjr = Ljr dt while respecting the divergence-free current constraint, it is 
easier to vary 

+ ■Jg))dt (3.311) 

Jti 

using unconstrained variations. Here the Lagrange multiplier ^ is a time-dependent function on Q 
modulo constant functions. 

The Euler-Lagrange equations associated with the action S'jr are given as follows. 


Varying gg — Varying the multi-species phase space fluid configuration gives 



+ iv;d^)^ -f egdiy ■ XG - v ■ At{J)) + gs* 


~ 0 - 


(3.312) 
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Here J = J{V, F) (note that J ^ Jol) and 5g^F is a 1-form on the single-particle phase space 
that can be thought of as the functional derivative of F with respect to Qs- To be precise, 


d 

de 


F{gse) 

0 




(3313) 


where = dgs o g^ In deriving the expression (3.312), it is useful to make use of the self-adjoint 
property of the transverse vector potential. 


{Ji,At{J2)) - (Ar(Ji), J2), 


(3.314) 


where Ji, J 2 are arbitrary vector fields on Q. It is also useful to note that 3.312 implies that the 
velocity Vs can be decomposed as 


K = ujTs + FsiesAjr ■ dx), 


(3.315) 


where is the Poisson tensor associated with the symplectic form Us = — and 


^Ts Ps(<^gs'T' -\- 6^d['U ■ CgA'j'i^Jjr^ gsif 




^s) 


(3.316) 


Varying X — Varying the scalar A gives 


6F 


(3.317) 


Note that 6F/6X lives in the dual to the space of scalar functions modulo constant functions, 
which is precisely the space of functions on Q with vanishing integral. 
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Varying Fos — Varying the reference phase space density Fos gives 


Xs + ■ {At{J) - Vi^)] - = 0. 


(3.318) 


This equation implies that the phase space function Xs behaves very much in the same way as the 
phase of a quantum wave in the WKB approximation. 

Varying Xs — Varying the phase-like function Xs gives 



(3.319) 


Varying — As expected, varying the Lagrange multiplier -0 gives 


V • Jo = 0. 


(3.320) 


Taken together, these Euler-Lagrange equations comprise a puzzle that must be solved in or¬ 
der to find the components of the Hamiltonian vector field Xjr. Some parts of the puzzle are 
simple. For instance, the time derivative of Fos is given immediately by Eq. (3.319). However, 
Eq. (3.312) represents a much more serious challenge. There the velocity field Vs is woven 
into various terms in a rather intricate manner. Complicating matters further is the fact that the 
Eagrange multiplier xjj must somehow be eliminated from the equations. 
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The key to solving this puzzle is obtaining expressions for the quantities 

= VA - At[J{Xo, LvF)] (3.321) 

Aot = Vij- At[J{V, F)], (3.322) 


in terms of functional derivatives of F. In principle, this task is not so difficult. However, obtaining 
the desired expressions in a useful form is a subtle enterprise. It turns out that a particularly useful 
way of writing the expressions is 


A^ = [n- 1] 
Aot =[n-l] 


[i-oJT + /ioVG 

a;2/c2 




Here the various current densities are given by 


(3.323) 

(3.324) 


J^ = J(F(6,F),F) 

(3.325) 

J^^ = J(iy^,F) 

(3.326) 

= J{X,,g,{5F/5x)) 

(3.327) 


and the /-dependent linear operators n, jj, are given by 


n = —V 


-1 


rAHi 




co„ 

/i = 1 - -icnr. 


(3.328) 


(3.329) 
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As we will see, the useful property of this representation is that the operator n{(?/oj^) is self- 
adjoint. The self-adjoint property ean be verified by first noting Pl = '^V^, whieh implies 

C CJp 

2 2 

Then one uses the fact that AIIt commutes with the operator V = 1 — Allr to conclude that 

LOp 

{n-^y = n^. It is also useful to be aware of the identity 

= -V-\ (3.331) 

With equations (3.323) and (3.324) in hand, we can now express the components of the Hamil¬ 
tonian vector field Xjr as 



Vs = vjrs + f‘s{esAjr-dx) (3.332) 

Fos = (3.333) 

OXs 

A ^ 

As = - giy&sivTs) - esV ■ Aojr) (3.334) 

x = g[v-A^]. 


(3.335) 





Chapter 3. Boundary terms and Poisson brackets 


90 


Upon substituting these expressions into the identity [Q, = Lxj,G, we find the following ex¬ 

pression for the Poisson bracket, 


le.-F] + ‘'SFJSx) 




oji 


Jr + VG[5J^/5\] ) + ( ^^oKJTx) 




Jg + VGidG/SX]^ 

( 3 . 336 ) 


This Poisson bracket is defined on a space larger than the phase space we are actually interested 
in. The “unimportant” variables are the Xs and the Qs- We can find the bracket on (F, A)-space 
using Poisson reduction. First we will reduce by the abelian group Yhs G°°{TQ), which acts by 
translation, Xs Xs + Sxs- H is simple to verify that both the bracket and Hamiltonian given 
earlier are invariant under this group action. The first reduced bracket is therefore 


[G, J^]ri = Fs), ( 3 . 337 ) 

5 

where F, G are functionals of {g, Fo, A) and 


+ esd[v ■ Aor])- ( 3 . 338 ) 

Next we will reduce by the product of diffeomorphism groups nsDiff(T(5) which acts (on the 
right) according to {gs, Fos) ((?« o hg, h*gFos). Again, the bracket and Hamiltonian given earlier 
are invariant under this group action. The second reduced bracket is 


[G,F]r2 = J2(i^F.+ CgV ■ Aog, Ffs + ^gV ■ Aor}, F), 


(3.339) 
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where 


Ao:f = F) + fx,VG[5J^/5\]^, 

and J^F = ST/dF. Here, F, Q are functionals on (F, A)-space. 


(3.340) 


3.6 Force-free electrodynamics 

Force-free electrodynamics is a single-fluid model sometimes used in astrophysics to describe 
extremely-magnetized relativistic plasmas. As discussed in Gralla and Jacobsen (2014), the defin¬ 
ing equations (which are Lorentz invariant) are given by 


dtB = —cV X E 
dfE = cV X B — AttJ 


J = 


AttB^ 

E-B = 0. 


(V ■ E)E xB + {B-VxB-E-Vx E)B 


(3.341) 

(3.342) 

(3.343) 

(3.344) 


3.6.1 Phase space variational principle for FFE 

For the sake of deriving the FFE Poisson bracket by reduction of a non-degenerate bracket, it is 
useful to formulate a phase space variational principle on an augmented phase space P+. The space 
P+ is the space of all ordered lists (P, A, g, po, x) where 

• P is the electric field 1-form 

• A is the vector potential 1-form 

• is a diffeomorphism of that represents the fluid configuration 
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• po is the reference charge density 3-form 

• X is an auxiliary scalar field, 
and E and A are constrained to satisfy 


EAdA = 0, 


(3.345) 


which guarantees that the electric and magnetic fields are perpendicular. The Lagrange 1-form on 
this space is given by 


Q{E,A,g,po,x) = {A®p,u/c) - ^{E,A/c) + (po,x)- (3.346) 


Here u = g o g angle brackets denote a “natural pairing”, p = g^Po, and parentheses denote the 
standard inner product of differential forms on The Hamiltonian functional is given by 


A, g, Po, x) = ^i.E, E) + ^{dA, dA). 
Stt Stt 


(3.347) 


The phase space Lagrangian is given by 


Ln = Q{E, A, g, po, x) - 'H{E, A, g, po, x)- 


(3.348) 
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To verify that this Lagrangian reproduces the equations of FEE, it is useful to calculate the first 
variation of Ljr, which is given by 

= ( » k/oPi + Tj3/c+ - sj .m) 

+ [-^/c - t^u/cB - {A/c)g^ {po/po)] ® P,^ 

+ {9*{A{u/c)) + x,5po) 

+ {-po,Sx) 

+ ^ e { 6 E , 6 A , 6 g , Spo , Sx ). (3.349) 

at 

Here B = dA, \B\ = ^/*{B A *B), b = *B/\B\, the dot product between two /c-forms a, (5 is 
given by 


a ■ (3 = *{a A */3), 


(3.350) 


and if a is a 1-form, a± = a — {a ■ b)b. Note that 6E is constrained to satisfy 
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i.e. the parallel variation of E is eompletely determined by the variation of B. By setting the 
variation of equal to zero, the following equations of motion emerge: 


Po = 0 (3.352) 

X = -g*{A{ulc)) (3.353) 

A/c = -E (3.354) 

6B = *[ iu / cP ] + (3.355) 

E = iu / cB , (3.356) 

whieh are equivalent to the equations of foree free eleetrodynamies provided An * p = — 5 EdX 
t = 0 (note that 47r * p + <5i? is a eonstant of motion). 


3.6.2 Symmetries of the phase space Lagrangian 

Symmetry under translations in x 

Consider the abelian group Go = C°°(M^), whose typieal element we will denote r. Go acts on 
the augmented FFE phase space according to x X + '?'• The phase space Lagrangian is invariant 
under this symmetry, which implies that 0(0, 0, 0, 0, r) = (po, r) is a constant of motion for each 
T G Go- Because r is an arbitrary smooth function, po must be independent of time. 
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3.6.3 Symmetry under time-independent gauge transformations 

Consider again the abelian group Go- Let s G Go be a typical element. This group acts on P+ in a 
second way, namely 


A + ds (3.357) 

X ^ X - - 9 * s - (3.358) 

c 

Note that this transformation preserves the constraint E A dA = 0. The phase space Lagrangian is 
invariant under this symmetry which implies that 

/is = 0(0, ds, 0, 0, —g*s/c) = -- —(dE + drr * p, s) (3.359) 

drrc 

is a constant of motion for each s E Go- Thus, 5E + drr * p is a constant of motion. 

Symmetry under particle relabling 

Let Gi = Diff(M^) be the non-abelian group of diffeomorphisms of with typical element 
h E Gi- There is a (right) Gi action on P+ given by 


goh 
Po H- h* po 

X h*x- 


(3.360) 

(3.361) 

(3.362) 


The phase space Lagrangian is invariant this symmetry. Thus, the contraction of the Lagrange 
1-form with the infinitesimal generator of this symmetry is a constant of motion. The infinitesimal 
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generator is given by 


= {5E, 6A, 6g, 6p, Sx) = (0,0, Tg o L^x), (3.363) 

where ^ is an arbitrary veetor field. This implies that 

At? = 0(^P+) = {{9*^/c + dx) O Po, 0 ) (3.364) 

is a eonstant of motion for eaeh Beeause ^ is arbitrary and po is a eonstant of motion, the quantity 
Ao = g*A + dcx is a eonstant of motion. To understand this eonservation law, consider the exterior 
derivative of Ao, 


dAo = g*B. (3.365) 

By applying the pushforward by g to each side of this expression, and noting that Ao is a constant 
of motion, we conclude that the magnetic field 2-form is advected by the fluid velocity u, 

B = B = —LuB. (3.366) 

But this implies that the vector potential A is advected modulo an exact 1-form, 

A = g^:Ao + d^jj. (3.367) 

Thus, a second conclusion we can draw about the conservation of is that dcx can be interpreted 
as the difference between A in the temporal gauge and A in the “advection gauge”, wherein A is 


advected as a 1-form. 
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Symmetry under isometries of configuration space 

Let G 2 be the non-abelian group of isometries of with typical element R. G 2 is naturally a 
matrix group, and so it is equipped with a natural (left) action on 


X I—)■ Rx. 


(3.368) 


Given an element z/ G 02 we can therefore define an infinitesimal generator on R^, 




de 


exp(ez/), 


(3.369) 


where exp denotes the matrix exponential. The left action of G 2 on R^ lifts to a left action on 
given by 


y4 1 —y 

(3.370) 

E R^E 

(3.371) 

g^ Rog. 

(3.372) 


The infinitesimal generator of this action is given by 

z/p+ = {6E, 5 A, 6g, 6p, 6x) = -L^^A, uq o g, 0,0). (3.373) 

The phase space Lagrangian is invariant under this symmetry. Therefore 

Pu = Q{^p+) = P + SE) + ^(P,6^qP) 


(3.374) 
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is a constant of motion for each u E Q 2 - This conservation law is equivalent to the global conser¬ 
vation of linear and angular momentum. Provided initial conditions are chosen so that the Gauss 
equation is satisfied, is equal to the total momentum of the electromagnetic field. 

Combining the symmetry groups Go, Go, and Gi. 

There is a “big group” that encodes three of the symmetries discussed so far: phase translation, 
gauge transformation, and particle relabling. The big group will be denoted G = Go x {Go x Gi), 
which is suggestive of the fact that H is a direct product of the abelian gauge group Go with the 
non-ablian semidirect product Go x Gi. The identity element is e = (0, 0, id) E Go x Go x Gi. 
The group product is given by 

(si, Ti, hi) * (S2, T2, h2) = (si + S2, h^Ti + T 2 , hi O h2), (3.375) 

which can be seen to be associative by direct calculation. The group inverse is given by 

(s, r, h )~^ = (—s, — h ^ T , h ~^). (3.376) 

While the structure of this group may appear a bit mysterious, it can be uncovered in a straight¬ 
forward way as follows. By an abuse of notation, denote the right actions of Go, Go, and Gi by 
Rs, Rt, and Rh, respectively. We have 


Rs{E, A, g, po, x) = {E,A + ds, g, po, x - 9 *s/c) (3.377) 

Rr{E, A, g, Po, x) = {E, A, g, po, X + u) (3.378) 

Rh{E, A, g, Po, x) = {E, A,goh, h* po, h*x)- (3.379) 
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Note that we have the following eommutation relations amongst these maps 

(3.380) 

(3.381) 

(3.382) 


RgO Rt = Rt ° Rs 
Rg o Rfi Rfi o Rg 

Rh ^ Rt Rh*T ^ Rh' 


Therefore we have the identity 

(yRs2 ® Rt2 ® ® (-^Si ^ Rt\ ® ^h\) Rs\+S2 ® Rh2T\+T2 ® Rhioh2 • (3.383) 

This identity strongly suggests defining the group produet given above. Morevoer, by eombining 
this identity with the faet that the produet (3.375) does indeed satisfy the group axioms, we ean 
quickly deduce the appropriate right action of ^ on P+, namely 

R{s,T,h) = RsO Rt° Rh- (3.384) 

The identity (3.383) can be rewritten as 

R(s2,T2,h2) ° R{si,Tl,hl) R{si,Tl,h-i)*(s2,T2,h2)l (3.385) 

which is the most non-trivial property satisfied by a right Q action. The other properties of a right 
action follow immediately from the definition (3.384) and the commutation relations. 

The quotient space P = P+/Q 

The Poisson bracket, {•, on P+ given by inverting the symplectic form —d© will automat¬ 
ically be symmetric under the group Q because 0 has the same symmetry. Therefore there is a 
reduced Poisson bracket on P = P+/Q that is given formally as follows. Let tt : P+ — )■ P 



Chapter 3. Boundary terms and Poisson brackets 


100 


be the projection map that sends a point in P+ to its orbit under the action of Q. Given a pair 
of functionals F, G : P —)■ M, we can pull them back to P+, thereby obtaining the functionals 
71*F, 7i*G : P+ —)■ M. Each of these functionals on P+ is automatically ^-invariant, which implies 
that the functional {7r*P, tt*G}p^ is also ^-invariant. Because tt is surjective, there is therefore a 
unique functional, {P, G}p : P —)■ M, that satisfies 

7r*{P, G}p = {ti*F, 'k*G}p^. (3.386) 

The formula (3.386) defines the Poisson bracket on P. The Jacobi and Leibniz identities are 
straightforward to check. 

Because P represents a “physical” phase space for FEE, the bracket on P is in some ways more 
desireable than the bracket on P+. It is therefore useful to have a concrete model of the space P. 
The purpose of this subsection is to supply this model. 

It will turn out that P can be represented as the space of triples {E,B,p), where E is the 
electric field 1-form, B is the (exact) magnetic field 2-form, p is the charge density 3-form, and the 
electric and magnetic fields are constrained to satisfy P A P = 0. To see this, first note that there 
is a surjective map tt : P+ —)■ P given by 

A g, po, x) = {E, dA, g^po). (3.387) 

Next suppose (P, A, g, po, x) and (P', A', g\ x') each map to (P, P, p) under tt. Immediately 
we see that E' = E and dA' = dA, which implies that there is some s & Go such that A' = 
A -|- ds. If we now set h = g~^ o g' and r = x' — h*x + h*g*s/c,a. simple calculation shows that 
{E', A' , g' , p'^,x') = R{s,T,h) (P, A, g, po, x)- This shows that points in tt ^(P, P, p) are all on the 
same ^-orbit. Because it is also true that tt o R[s,T,h) = tt for each (s, r, h) G Q, we must therefore 
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conclude that the preimages of points in P along vr are preeisely the ^-orbits. This proves that 
P+/Q is diffeomorphic to P. 

3.6.4 Derivation of the Poisson bracket on P 

In order to eompute the Poisson braeket on P, we will eompute 71*0} p^, where F,G : P ^ 

M are funetionals on P. By the definition of a Hamiltonian veetor field, we have 


{7r*P,7r*G}p^ = (d7r*P)(X,*G). 


(3.388) 


In order to evaluate the right hand side of Eq. (3.388), we must eompute the eomponents of the 


Hamiltonian veetor field Xt,*g = {Eg, Ag, Qg, PoG, Xg)- This eomputation is most easily done by 
analyzing the Euler-Lagrange equations assoeiated with the phase spaee Eagrangian 



The first variation of Lg is given by 



+ ( {-Ajc - iujcE - {Alc)g^{polpo) - d5G/5p] ® p,^ 


{h-[A/cPATi5G/5E])E 

\B\ 



+ { 9 *{A{u/c)) + x- g*6G/6p, 6po) 


+ {-po,Sx) 

+ ^9, Spo, Sx)- 


(3.390) 
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It follows that the components of the Hamiltonian vector field Xt,*g are given by 


Po = 0 

-Ag = —dvrcKc-x — 47rd-— 
c op 

1 ■ 47r*p , f E J *E 5G\\ 

-JG = -^*{oLGAb)- \^—-d{*cxG±) + b-5\^aG\\j^^-— j jb 

Ip. 47r*p . J *E SG\ (E 

=-* {Attolg f\b)- Ati5 - — J + ■ S{*ocg±) j b 

Xg = 9 *(^^ - 41 ( mg / c )^, 


(3.391) 

(3.392) 

(3.393) 

(3.394) 

(3.395) 


where cxg = SG/6E — {l/4:7i)d6G/6p and jg = *GgP is current density 1-form. 

After some tedious but straightforward calculations involving substituting the previous expres¬ 
sions into the identity (3.388), we find that the Poisson bracket on P is given by 


{E 


, . f {4:71 * p)B \ 

, G}p = -dvrcf —1^1^—, cxp A aGj 

+ 47Tc(aF±,s(^ - ttGll *e ) ) -47rc(Q:G±,5( ^ - 


CtF\\ * 6 


(3.396) 


where e = E/\B\. 

It is not difficult to show that any functional of dE -f dvr * p is a Casimir of the bracket {•, -jp. 
Therefore the submanifold of P defined by 5E -f 47r * p = 0 is a Poisson submanifold that can 
be parameterized by the space Po of pairs {E^ B) that satisfy i? A i? = 0. Being a Poisson 
submanifold, Po has a bracket {•, -jpo that is naturally induced by {•, -jp. The expression for this 
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bracket on {E, 5)-space is given by 


{F,G}p^=4nc 


-|-47rc 


—47rc 


( 


{dE)B 

6E 

SG\ 

\B? ’ 

Ie ^ 

se) 



SG 

SE±^^^ 

[Jb~ 

' ^11 


f6E 

SE 

SE± ' 

[Ib 

'SE\\ 


* e 


* e 


(3.397) 




Chapter 4 

Energetically-consistent gyrokinetic 
collision operator 

4.1 Introduction 

One of the greatest unsolved problems in the theory of magnetieally-eonfined plasmas is un¬ 
derstanding and eontrolling the turbulent flux of partieles and heat into a fusion reaetor’s 
wall Kikuehi and Azumi (2012). It is believed that the predominant eause of these fluxes is 
low-frequeney fluetuating eleetromagnetie fields with wavelengths on the order of the gyroradius. 
While a eollisionless gyrokinetic model of these fluctuating fields has been developed that is fully 
consistent with the First Law of Thermodynamics (for a recent review see Brizard and Hahm 
(2007)), this energetically-consistent model has the serious flaw of ignoring collisions altogether. 

In order to accurately describe irreversible plasma transport processes, the effects of colli¬ 
sions must be incorporated into gyrokinetic theory. Previous work on linear gyrokinetic collision 
operators Abel et al. (2008); Li and Ernst (201 1); Madsen (2013a) assumed a strict two-scale sepa¬ 
ration between a large-scale equilibrium distribution function Fo and a small-scale fluctuating part 
5F = F — Fo. Conservation properties of the collision operator in Abel et al. (2008), for example. 
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were discussed in the gyroBohm limit. Here, we will focus on nonlinear gyrokinetic collision op¬ 
erators for a global full-F approach that do not make this split, and that can thus investigate more 
completely the possible effects of finite e = pi/L in experiments, such as corrections to gyroBohm 
scaling and non-local turbulence spreading (see footnote 5 on p. 427 in Brizard and Hahm (2007).) 

When finite-e effects are accounted for, preserving exact conservation properties, and therefore 
ensuring consistency with the First Law of Thermodynamics, is a nontrivial unsolved problem. 
The collision operators in Abel et al. (2008); Li and Ernst (201 1), for example, were obtained by 
transforming a particle-space collision operator with exact conservation properties into the lowest- 
order guiding center coordinates. While this approach guarantees the existence of energy and 
momentum-like quantities that annihilate the collision operator, these same quantities are not con¬ 
served by the full-F collisionless gyrokinetic system, and therefore fail to be conserved by the 
full-F collisional system. More generally, existing gyrokinetic collision operators are not ener¬ 
getically consistent in a full-F formalism because: (a) the gyrocenter coordinate transformation, 
and therefore any collision operator transformed into gyrocenter coordinates, is only known as an 
asymptotic expansion in the gyrokinetic ordering parameter e; and (b) replacing the asymptotic 
expansion of such an operator with a truncated power series destroys exact conservation laws. The 
purpose of this Chapter is to present the first collisional formulation of global full-F gyrokinetics 
with exact conservation laws. 


4.2 Electrostatic Model 

For the sake of simplicity, our discussion will focus on quasi-neutral electrostatic gyrokinetics (for 
instance, see PaiTa and Calvo (201 1)). However, the ideas behind our discussion apply equally- 
well to electromagnetic gyrokinetics (for example, see Sugama (2000).) Our primary result con¬ 
sists of an expression for the non-linear Landau operator in gyrocenter coordinates that is cor¬ 
rected by small terms to ensure exact energy and momentum conservation [see Eq. (4.24).] These 
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correction terms are analogous to the -denominators in the Hamiltonian guiding center theory 
introduced by Littlejohn Littlejohn (1981); they do not increase the theory’s order of accuracy, but 
they are essential to include for the sake of ensuring exact energy and momentum conservation. 

As a first step, we review how the energy conservation law is discussed in collisionless ki netic 
theory. The governing equations of collisionless electrostatic kinetic theory are the Vlasov-Poisson 
equations. 


dtf, + {fs,Hs}=0 (4.1) 

Aip = -47rp(/), (4.2) 

where fs is the species-s distribution function, (p is the electrostatic potential, p{f) is the charge 
density, Hg = p‘^/2ms + CsP, and {•, •} is the standard canonical Poisson bracket. Equations 
(4.1)-(4.2) conserve the total energy 

where {■,■) denotes the standard -pairing of functions on configuration space and dz = dx dp. 
Because binary collisions conserve energy, Eq. (4.3) must also be conserved in collisional kinetic 
theory. In particular, if the Vlasov-Poisson equations are modified by the addition of a bilinear 
collision operator. 


dtfs + {fs.Hg}=Y,Cs-sUsJs) 
s 


(4.4) 


Ap = -47rp(/), 


(4.5) 
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then Css niust be ehosen to satisfy the eondition 

s ' ' 

= J2[ HsCs,{fs,h)dz. (4.6) 

S,S 

Because this identity must hold for an arbitrary multi-species distribution function, the collision 
operator therefore has to satisfy the well-known identities 

j HsCs-sifs, /.) dz + j H-sC-ssif-s, fs) dz = 0 , (4.7) 

which express the fact that the energy gained by species s due to collisions with species s is 
precisely the energy lost by species s due to collisions with species s. The non-linear Landau 
operator (summation rule is implied), 


= (4,8) 

satisfies the identities (4.7), and therefore defines an energetically-consistent collisional kinetic 
theory. Here In A; the 3-component vector 7 ^^ is 

= j d{x- x) ^"\z, z) Ass{z, z) dz] (4.9) 


the 3x3 matrix is given by 




(4.10) 
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where P(^)=I — the orthogonal projeetion onto the plane perpendieular to the veetor the 
veloeity difference is given by 

Ws,{z,z) = {x,H,}{z) - {x,H,}{z)- (4.11) 


and the vector 


Ass{z, z) = fs(z){x, f^}(z) - {x, fs}{,z)f-s{z). (4.12) 

When comparing this form of the Landau operator to more conventional expressions, it is useful 
to note that {x, g} = dpQ, where g is any function on phase space, so that the collision operator 
(4.8) describes collisions in momentum space. Moreover, the identities (4.7) follow immediately 
from the fact that the velocity difference is a null-eigenvector of the matrix 

4.3 Electrostatic Gyrokinetic Model 

In order to apply this same argument to gyrokinetic theory, we start with the gyrokinetic Vlasov- 
Poisson system 


dtF, + {F,,Hf}f = 0 (4.13) 

V-P = p{F). (4.14) 

Here, Fg is the gyrocenter distribution function; g) is the electrostatic potential; {•, is the guid¬ 
ing center Poisson bracket; 

Hf = Hf + e. W + i}r> = K,(E) + 


(4.15) 
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is the gyrocenter Hamiltonian; •^( 2 ;) = ip{X + pos), where pos is the lowest-order guiding-eenter 
gyroradius; {■) denotes the gyroaverage; ^ denotes the gyroangle antiderivative of ip = 7p — {ip); 
Ks{E) is the gyrocenter kinetic energy; P = — 51C/5E h the gyrocenter polarization den¬ 
sity; ^ si PsKs{E) dzf; and dzf denotes the guiding center Liouville volume element. 
These equations govern collisionless quasineutral electrostatic gyrokinetic theory in the “high- 
flow” regime (see Krommes and Hammett (2013) and references therein) and they conserve the 
total energy, 



(4.16) 


exactly. Note that the quasineutrality equation (4.14) implies that this system governs plasma dy¬ 
namics on time scales long compared to the period of plasma oscillations. 

The equations governing collisional gyrokinetic theory are given by adding a bilinear collision 
operator to the gyrokinetic Vlasov-Poisson equations, 



(4.17) 


S 


V ■P = p{F). 


(4.18) 


Because the conservation laws of ordinary collisional kinetic theory are consistent with those of 
collisionless kinetic theory, the gyrokinetic collision operator must not alter the conservation 
of^^T Thus, 



(4.19) 
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This identity will be satisfied for a general multi-species gyrocenter distribution function if and 
only if 



(4.20) 


which is the gyrokinetic version of Eq. (4.7). The identities (4.20) must be satisfied exactly by any 


energetically-consistent gyrokinetic collision operator. 

4.4 An energetically-consistent collision operator 

While Eq. (4.20) imposes important qualitative constraints, they cannot determine the form of the 
gyrokinetic collision operator by themselves. A quantitative constraint is necessary as well. To 
this end, it is important that the gyrokinetic collision operator agrees with the the transformation of 
the particle-space Eandau operator ^ into gyrocenter coordinates, at least up to some desired order 
in the gyrokinetic ordering parameter e. Is it possible to satisfy these qualitative and quantitative 
constraints simultaneously? The answer is “yes”. 

We have discovered an accurate gyrokinetic collision operator that is consistent with the con¬ 
servation laws of collisionless gyrokinetic theory, and therefore the first law of thermodynamics. 
The form of the operator is suggested by the somewhat-peculiar presentation of the particle-space 
Eandau operator given earlier. Eet ys = X + pos and define the gyrocenter velocity difference 



(4.21) 


'Necessary conditions for the use of the Landau operator are uic < ojp and {dtF)/{u!pF) < 1. When these 
conditions are not satisfied, our discussion must be modified. 
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the associated 3x3 matrix 


“gy 






-P[EEf(z,x)], 


(4.22) 


and the vector 


AlUz,z) = F,(z){y„F,}r(z) - {y., Fjf (x)F,(x). 


(4.23) 


The energetically-consistent gyrokinetic Landau operator is given by 




(4.24) 


where 


7 gy(z) = / 6l^(z,z)Q"y(z,z)All(z,z)dzf, 


(4.25) 


and ()f|(x, z) = d(ys(z) — ys(z))- Note that this operator depends explicitly on the electric field 
through the gyrocenter Hamiltonians that appear in Eq. (4.21). Using a straightforward, but tedious 
argument that is not reproduced here, we have shown that this operator agrees with the Landau 
operator transformed into gyrocenter coordinates with leading-order accuracy. 

Because the proof is simple, we will now show explicitly that the gyrokinetic Landau-Poisson 
system (4.17) defined in terms of the collision operator (4.24) has exact conservation laws for en¬ 
ergy and momentum. We hope to convey the similarity of this demonstration with the analogous 
demonstration for the ordinary Landau-Poisson system (4.4)-(4.5). However, a word of caution 
is in order here. It is essential that the guiding center Poisson brackets that appear in Eq. (4.24) 
be genuine Poisson brackets (i.e., the brackets must satisfy the Leibniz and Jacobi identities). 
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Dropping terms from a bracket that satisfies these properties will destroy the gyrokinetic Landau- 
Poisson system’s exact conservation laws. 

4.5 Energy conservation 

Proving that the gyrokinetic Landau operator (4.24) satisfies the identities (4.20) is very similar 
to proving that the particle-space Landau operator satisfies the identities (4.7). Setting Sgs = 
f Fg) dzf, it is simple to verify that 

+ = (LLf )tQ--Ai^f| dzf dzf, (4.26) 

where all two-point quantities in the integrand are evaluated at {z, z) and denotes the ordinary 
matrix transpose. Because Q|y is a symmetric matrix with null eigenvector Wff the right-hand- 
side of this equation vanishes exactly. Thus the gyrokinetic Landau operator (4.24) satisfies the 
identities (4.20) exactly, and the gyrokinetic Landau-Poisson system (4.17) has an exact energy 
conservation law, dS^^/dt = 0. 

4.6 Toroidal momentum conservation 

We will prove that if the background magnetic field is axisymmetric, then the gyrokinetic Landau- 
Poisson system conserves the total toroidal momentum 


= P^sFsdzf, 


(4.27) 
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where p^s is the guiding eenter eanonieal toroidal momentum If the baekground magnetie field 


has additional symmetries, a similar proof of the eonservation of the eorresponding total momen¬ 
tum can easily be constructed. The time derivative of Eq. (4.27) yields 



(4.28) 


where is conserved exactly by the gyrokinetic Vlasov-Poisson system. Here, we find 




(4.29) 


Now using the fact that p^s is the generator of infinitesimal toroidal rotations, we can see that 
{Vs ^P 4 >s}f = X Vs, where is the unit vector along the axis of rotation. Therefore the vector 
quantity {{Vs^P^sjf - {ys,P<t>s}f) (iS = x {ys - ys) = 0, which follows from standard 
^-function properties. This shows that P^ss + P<i>ss = 0> which in turn implies total toroidal 
momentum conservation dP^/dt = 0. 

4.7 Entropy production 

As we have discussed, these conservation laws ensure that the gyrokinetic Eandau-Poisson system 
is consistent with the the Eirst Eaw of Thermodynamics. On the other hand, they do not directly 
imply that the gyrokinetic Eandau-Poisson system is consistent with the Second Eaw of Thermo¬ 
dynamics. To verify that entropy is indeed a non-decreasing function of time, we have computed 

^Rather than give an explicit expression for , which will depend on ones choice of guiding center representation, 
it is better to define it operationally via the guiding center Poisson bracket; for each phase space function /, the 
canonical toroidal momentum satisfies {f,P(ps}f = d^f, where is the toroidal angle derivative. 
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the time derivative of S' = —J2s I dzf and found 


dt 2 


// 



(4.30) 


Beeause Qgy is a positive semi-definite matrix and the distribution funetion is positive the right- 
side of Eq. (4.30) is non-negative, which is the desired result. 

Note that this proves one “half’ of a gyrokinetic version of Boltzmann’s if-theorem. The 
missing ingredient is a complete characterization of the distributions that satisfy dS/dt = 0, i.e. 
the gyrokinetic Maxwellians. Because the guiding center Poisson bracket is rather complicated, 
we have not yet found a complete characterization. However, we have verified that the distribution 



(4.31) 


where Zg = J exp(—if|^/T) dzf is the partition function, maximizes the entropy. We leave the 
characterization of the most general gyrokinetic Maxwellian, which would be useful for the sake 
of deriving dissipative gyrofluid models with exact conservation laws Madsen (2013b), as a topic 
for future study. 

4.8 Gyroaveraging 

When the collision frequency is much smaller than the gyrofrequency Brizard (2004), the full 
gyrokinetic Eandau operator (4.24) can be replaced with that operator’s gyroaverage, (Cf|). When 
this is done, the gyrokinetic Eandau-Poisson system becomes the gyroaveraged Eandau-Poisson 

^Positivity of the distribution function is also guaranteed by the positive semi-definiteness of Q|y. 
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system, 


dtF, + {F„ Hf’}f = '^{CSiFs, Fi)) (432) 

5 

V-P = p{F), (4.33) 

where Fg is now interpreted as the gyroaveraged part of the distribution function. Because the 
functions Flf and p^s are independent of the gyrophase, the proofs of energy and momentum 
conservation given earlier work with C% replaced by (Cf|). Thus, the gyroaveraged Landau- 
Poisson system has exact energy and momentum conservation laws. 

4.9 Linearization 

Closely related to the gyroaveraged Landau-Poisson system is the collisionally-linear gyroaver¬ 
aged Landau-Poisson system. 


a,F. + {F., }f = ^ (iC“ + icr), 

S 


(4.34) 


V-P = p(F), 


(4.35) 


where the linearized test-particle and field-particle collision operators are 

6Cfi\Fg) = (Cf|(F„ Fms-)), (4.36) 

6C^f{F,) = {CIUFms, F,)). (4.37) 


This system of equations is obtained from the gyroaveraged Landau-Poisson system by assuming 
Fg = Fms + 5Fg and then dropping the non-linear term in the collision operator, {Cll{5Fs^ SFg)). 
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Note that {CsI{Fms, Fms)) = 0 Because the gyrokinetic Landau operator satisfies the identi¬ 
ties (4.20), it is straightforward to prove that these equations have the same conservation laws for 
energy and momentum as the gyroaveraged Landau-Poisson system. 

4.10 Concluding remarks 

The key to deriving an energetically-consistent formulation of collisional gyrokinetics was first 
expressing the particle-space Landau operator in terms of Poisson brackets “as much as possible,” 
which was an idea first championed by Brizard in Brizard (2004). In particular, the identity 

V - V = {x,Hs}{z) - {x,Hs}{z) (4.38) 

suggests that the appropriate definition of the gyrocenter velocity difference is given by Eq. (4.21). 
This idea, together with the procedure given earlier for determining the energetic consistency con¬ 
straints, appears to be appropriate for deriving energetically-consistent collision operators for other 
reduced plasma models as well. In future work, we will report on the energy-conserving collisional 
formulations of electromagnetic gyrokinetics and oscillation center theory. 

We note that, although the gyrokinetic Landau operator (4.24) and its linearized forms (4.36)- 
(4.37) may prove difficult to implement numerically, they identify the proper formalism for the 
inclusion of collisional transport in gyrokinetic theory. Hence, these gyrokinetic collision operators 
form the basis from which approximations can be implemented for practical applications. 

Lastly, by setting 99 = 0 in the above formulas, our results reduce to an energy-momentum- 
conserving guiding center collision operator. This operator would be ideally suited to incorporating 
collisions into orbit-following codes such as ORBIT White and Chance (1984); see Hirvijoki et al. 
(2013) for recent work on the Monte Carlo implementation of a 5D guiding center Lokker- 

"^Note that this identity does not contradict the message presented in Madsen (20 1 3a). In that reference, the gyroki¬ 
netic Maxwellian is dehned using only the lowest-order gyrocenter Hamiltonian. 
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Planck collision operator. All previous guiding eenter eollision operators that have been applied 
in orbit-following eodes either resort to ad hoc methods to ensure exact conservation laws 
Boozer and Kuo-Petravie (1981), or else do not fully aeeount for inhomogeneities in the magnetic 
field Tessarotto et al. (1994). 

4.11 Acknowledgements 

The results presented in this chapter were obtained in collaboration with Professor Alain Brizard. 
They ean also be found on the arXiv at arXiv:1503.07185. 



Chapter 5 

Hamiltonian mechanics of stochastic 
acceleration 

5.1 General theory 

5.1.1 Introduction 

The term “stoehastie aeeeleration” refers to the ehaotie motion of partieles subjeeted to a preseribed 
random force. Such motion occurs in myriad contexts; the turbulent electromagnetic fields present 
in the interstellar medium and the RF wave fields found in magnetic fusion devices are just two 
examples. In the astrophysical context, it is thought to be partially responsible for the presence of 
cosmic rays in our solar system Fermi (1949). In the magnetic fusion context, it might explain the 
presence of certain high-energy tails observed in the National Spherical Torus Experiment when 
neutral beams are fired into RF-heated plasmas Liu et al. (2009). 

Robust modeling of stochastic acceleration requires statistical approaches. The dominant 
approach is to employ the Fokker-Planck equation Sturrock (1966); Hall and Stun'ock (1967); 
Barbosa (1979); Petrosian and Liu (2004); Hamilton and Petrosian (1992) for the one-particle 
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distribution function. However, when studying Riehardson dispersion Riehardson (1926); 
Jullien et al. (1999), and more generally any phenomenon governed by the two-partiele distri¬ 
bution function Mukhopadhyay and Goree (2012), the one-particle Fokker-Planck equation is 
insuffieient. This is beeause spatial eorrelations in the random foree field prevent the two-partiele 
distribution funetion from faetoring as a produet of one-partiele distribution funetions. A superior 
statistieal model when multi-particle statisties are in question would be a Langevin equation for 
partiele trajeetories. A wisely-ehosen Langevin equation eould eapture the physies of the one- and 
two-particle distribution funetions while providing an attraetive means to perform Monte Carlo 
simulations of stoehastie aeeeleration. Currently, there are no satisfactory methods for finding 
sueh a Langevin equation. 

The purpose of this Chapter is to deseribe, for the first time, a systematie proeedure for pass¬ 
ing from a mieroscopie deseription of stochastic acceleration in terms of Hamiltonian equations 
of motion to the physically-correct Langevin equation for partiele trajectories in the long-time 
limit. We will also show that, aside from reproducing the correct multi-particle statistics, this 
Langevin equation inherits the Hamiltonian strueture of the mieroseopie dynamies. Speeifieally, 
we will show that the Langevin equation is a Hamiltonian stoehastie differential equation (SDE) 
Lazaro-Camf and Ortega (2008). Thus, this work proves that symmetries of the macroseopie phys¬ 
ical laws governing stoehastie aeeeleration lead to eonservation laws. 

We will foeus our attention on stoehastie aeeeleration problems similar to those studied in 
Sturroek (1966); Hall and Sturroek (1967); Barbosa (1979); Petrosian and Liu (2004). These eon- 
sist of a eolleetion of non-interaeting partieles moving through a preseribed Hamiltonian foree 
field. By assumption, the foree will eonsist of a small-amplitude perturbation superimposed over 
a time-independent baekground. The perturbed force felt by a particle will be assumed to have a 
correlation time mueh shorter than any bounce time associated with the perturbation, zero mean, 
and temporally homogeneous statisties. These assumptions preelude treating Coulomb eollisions 
beeause the polarization field produeed by a partiele eannot be modeled as a preseribed field; the 
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polarization force depends on the history of a partiele’s orbit. They also preelude the treatment of 
strong turbulence DuBois and Espedal (1978). 

5.1.2 The main idea 

Mathematically, this type of problem ean be described as follows. Eaeh particle moves through 
a 2n-dimensional single-particle phase spaee M aeeording to a dynamical law given by a time- 
dependent veetor field Xt; if Zt ^ M denotes the trajeetory of a partiele in M, then 


Zt = Xt{zt). (5.1) 

Beeause the only forees present are Hamiltonian, Xt must be Hamiltonian in the sense that there 
is some Poisson braeket {•, •} and some time-dependent Hamiltonian, Ht, sueh that i* = {z\ Ht}, 
where z^ denotes an arbitrary eoordinate system on M Grebogi et al. (1979). By standard mathe- 
matieal convention, this is written Xf = Xh^ Abraham and Marsden (2008). The presumed form 
of the force then implies Ht = Hq + eht, where e -C 1, i7o describes the mean time-independent 
baekground, and ht deseribes the small-amplitude random perturbation. Moreover, Xh^ evaluated 
on a partiele trajeetory must have a eorrelation time mueh shorter than some eonstant r, whieh, 
in turn, is mueh shorter than any bounee time assoeiated with the perturbation r^, Tac r Tb. 

Our goal in this Chapter is to find the correct coarse-grained version of the mieroscopic equa¬ 
tions of motion, Xnt- Speeifically, we seek a Langevin equation in the form 

fct = Xoizt) dt + J2 Mzt) 6Wt^ (5.2) 

k>l 

whose solutions correctly reproduce the late-time statistical behavior of solutions to the miero- 
seopie equations of motion. Here X^ are veetor fields on M that must be determined, are in¬ 
dependent ordinary Wiener proeesses, and 5 denotes the Stratonovieh differential Gardiner (2009) 
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(sometimes also written od). We will identify the Xk by demanding that Eq. (5.2) possess two 
properties: it must generate the Fokker-Planck equations for the one- and two-particle distribution 
functions, ft{z) andgt{zi, Z 2 ). The two-particle distribution function is defined such that the proba¬ 
bility particle 1 is in the region Ui <Z M and particle 2 is in the region f /2 C M at time t is given by 
fui fu 2 *^^ 2 , where d^; denotes the Liouville measure Abraham and Marsden (2008). Baxen- 
dale Baxendale (1984) has proven that a Langevin equation is uniquely determined by its one- and 
two- particle Fokker-Planck equations. Therefore, these conditions uniquely specify the Fangevin 
equation we seek. In particular, the requirement that two-particle statistics be accurately repro¬ 
duced is critical; Baxendale’s work implies that constraining the Fangevin equation only to be 
consistent with the one-particle Fokker-Planck equation would not identify it uniquely. 

Physically, the reason that the two-particle Fokker-Planck equation contains more information 
than the one-particle Fo kk er-Planck equation can be understood as follows. After a short amount 
of time At, the displacement of a particle initially located at zi at time t is given approximately 
by AtXt(zi). Similarly, the displacement of a particle initially located at Z 2 is nearly AtXt(z 2 ). 
Because the random force field generally has spatial correlations, Xt(zi) and Xt(z 2 ) are not statis¬ 
tically independent. Thus, the probability distribution of {z[, Z 2 ), where z' ^ zt + AtXt(zi), will 
not be given by the product of the distribution of with that of This failure-to-factor precludes 
determining the two-particle distribution function from the mere knowledge of the one-particle 
distribution function. Note that this is true in spite of the fact that these particles do not interact, 
because the random force is assumed to be prescribed, the time-evolution of zi is decoupled from 
the time-evolution of Z 2 . 
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5.1.3 Identifying the Langevin equation 

The one-particle Fokker-Planck equation associated with Eq. (5.2) is given by Gardiner (2009); 
Baxendale (1984) 


dt 


-div(/tXo) + ^ ^div(div(/iXfc)Xfc) 
^ k>l 


= AJt, 


(5.3) 


while the two-particle Fokker-Planck equation Baxendale (1984); Schmalfuss (2001); Kunita 
(1987) is given by 


.( 1 ) , .( 2 ) 

-2 ^div^^^div^^^ : gtXk{zi) O Xk{z2). (5.4) 

k>l 

The divergence operators in these expressions are defined relative to the Liouville volume form and 
the colon indicates the full contraction of second-rank tensors, a : b = a^^bij. Because these equa¬ 
tions follow from Eq. (5.2) via rigorous mathematics, we will refer to them as the mathematical 
Fo kk er-Planck equations. 

On the other hand, under our assumption that the correlation time of the perturbed force is 
much shorter than a bounce time, standard coarse-graining procedures Risken (1996); Bazant 
(2006) together with a decomposition theorem for time-ordered exponentials Earn (1998) lead 
to the late-time evolution laws for the one- and two-particle distribution functions associated with 
the microscopic equations of motion, Eq. (5.1). The physical one-particle Fo kk er-Planck equation 
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is given by 

Hq + — E[s2] I + —E[{{/i, Si} , Si}] 

= (5.5) 

while the physical two-particle Fokker-Planck equation (see the supplementary material for a 
derivation) is given by 

^ -f A^^^gt -f —E[a : (5.6) 

Ot T 

The notation introduced in these two equations is defined as follows: E denotes an expectation 
value; the functions si, S 2 are defined by 

51 = / exp{XXHo)*hr-xdX (5.7a) 

Jo 

52 = / / {exp(6X/iJ*h^_b,exp(aXHo)*/ir-a}d6da; (5.7b) 

^JoJo 

exp(F) M ^ M denotes the time-one advance map of the dynamical system defined by the 
vector field Y ; (exp(F)*h)( 2 ;) = h(exp(—F)( 2 ;)); the superscripts indicate which argument of gt 
that Ai and the exterior derivative d should be applied to; and a{zi,Z2) = E[X^^(2:i) (8)X5^(2:2)] 
is the two-point covariance tensor. 

The Xfc must be chosen so that the mathematical Fokker-Planck equations, Eqs. (5.3) and (5.4), 
are equivalent to the physical Fokker-Planck equations, Eqs. (5.5) and (5.6). However, a direct 
comparison of these two pairs of equations is difficult with Eqs. (5.5) and (5.6) in their current 
form. To eliminate this issue, we will obtain a special decomposition of the two-point covariance 
tensor a{zi,Z 2 ). 
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As a first step, notice that if we fix a one-form ^ e T*^M, then we can define a vector field 
on M by contracting ^ with a on the left according to 


Y^{z2) = 

= E[e(X,,(;^i))X,,(z2)]. (5.8) 

By forming all possible linear combinations of vector fields of this form, we can construct a (po¬ 
tentially infinite dimensional) linear space of vector fields Aronszajn (1950); Baxendale (1976), 
which we will denote "H, 


H = {linear combinations of ^ G T*M}. (5.9) 

Because each is of the form Y^{z) = X^iz) with H{z) = E[^(Xsi(zo))-Si( 2 :)], and the sum of 
Hamiltonian vector fields is again Hamiltonian, "H consists entirely of Hamiltonian vector fields. 
Moreover, following Baxendale Baxendale (1984, 1976), we see that "H is a real Hilbert space 
whose inner product is defined by the formula 

= m{Xs,{z^))'n{Xs,{z2))l (5.10) 

where ^ G T*^M and r] G T*^M. Therefore we may choose an orthonormal basis {ek}k>i for 
"H, where each must be of the form = X^^- A simple calculation then leads to the desired 
decomposition of a: 


a{zi,Z2) = '^Xh,{zi)(^Xh,{z2). 

k>l 


(5.11) 
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Using this decomposition of the two-point covariance tensor, it is straightforward to manipulate 
Eqs. (5.5) and (5.6) into the same form as Eqs. (5.3) and (5.4). After doing so, it is trivial to identify 
the correct Xk- Indeed, we have found that the physical Eangevin equation is given by 

Sz, = d( + A'sJj,) SWt, (5.12) 

k>l 

where 


Ho = Ho + -E[s 2], Hk = ^Hk (5.13) 

r 

Recall that the Xh^, are defined to be an orthonormal basis of the Hilbert space "H defined in 
Eq. (5.9). Also recall that all of the above manipulations have been performed under the assumption 
that the correlation time of the perturbed force felt by a particle is much shorter than any bounce 
time associated with the perturbation. 

Because the coefficients in the Eangevin equation for stochastic acceleration, Eq. (5.12), are 
all Hamiltonian vector fields, this equation is an example of a stochastic Hamiltonian system, the 
foundations of which are developed in Eazaro-Camf and Ortega (2008). It is in this sense that the 
Eangevin equation for stochastic acceleration inherits the Hamiltonian structure of the microscopic 
equations. In particular, SDEs of this type are known to arise from a stochastic variational principle 
for which Noether’s theorem applies. Thus, even at the dissipative macroscopic level, symmetries 
imply the presence of conservation laws. 

5.1.4 Example 1 

We will find the physical Eangevin equation for two example stochastic acceleration problems. 
Generally speaking, finding the coefficients of the physical Eangevin equation involves finding an 
orthonormal basis for the space "H, a task which may be analytically intractable. But, by Mercer’s 
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theorem Mereer (1909), this task ean be east as an eigenvalue problem for whieh there are existing 
numerieal solution methods. In any case, in these examples, the analytical route is tractable. 

First, consider a single-species, unmagnetized plasma subjected to a random weak electrostatic 
pulse at r-second intervals. Assume that the pulses are uniform in space and constant in magnitude, 
but uniformly and independently distributed in direction. Thus, the /c’th pulse is generated by a 
potential of the form (l)k{x,t) = {zk ■ x)(f)ou{t — kr), where Zk is a random vector uniformly 
distributed over the unit sphere and u{t) is a temporal windowing function localized att = r/ 2 . 

In order to find the Langevin equation governing the plasma dynamics at times much longer 
than r, we must (a) calculate si and S 2 using Eqs. (5.7a) and (5.7b), (b) find an orthonormal basis 
{XHi,}k>i for the space V. defined in Eq. (5.9), and (c) write down Eq. (5.12) with Hq and 
calculated using Eq. (5.13). The results of these three steps are as follows. 

(a) A quick calculation shows that 


(5.14a) 


Si = irioZ ■ X — miZ ■ v 


(5.14b) 


S 2 = const 


where mo = {q/m)(j)o u{s)ds, mi = (g/m) 0 o/ o^(t — s)M(s)ds, and g/m is the charge-to-mass 
ratio. 

(b) Each must be of the form = Xg^_^, where 


gi3-i{x, v) = -(mi/3 + rrioi) ' ("uri - rriox), 


(5.15) 


and f3 ,7 are arbitrary constant 3-component vectors. Using this expression, it is simple to find an 
orthonormal basis for "H. One is given by {Xjg^}fc=i.. 3 , with 



(5.16) 
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where {ei}i=i ,,3 is the standard basis for 

(e) Finally, the physieal Langevin equation is given by 


6x^ = n* dt H —^=mi (ifF* 
\/3r 

6v^ = 

vSr 


(5.17a) 

(5.17b) 


where z = 1, 2, 3. 

As is readily verified, the one-partiele Fokker-Planek equation for this SDE is given by 


^ + vV ft = ^{mlV'^ft + rriomiV ■ V^ft 
at or 


+ momiVy-Vft + mlVlft). (5.18) 


On the other hand, given an arbitrary funetion ((){x, v), the SDE 


6x^=v^dt+ ^ (cos(0) sin(0) 

vSr 

(5.19a) 

771 

5F- ° (cos(0) sin(0) , 

vSr 

(5.19b) 


where the W^’\ are six independent ordinary Wiener proeesses, will also generate Eq. (5.18). 
However, when 0 is not eonstant, the two-partiele Fokker-Planek equation generated by Eq. (5.19) 
will differ from the two-point Fokker-Planek equation generated by Eq. (5.17). This ean be verified 
using Eq. (5.6). The proeedure identified here seleets 0 = 0 as the physieal ehoiee. In partieular, 
it shows that a Eangevin equation with the eorreet one-particle Eokker-Planck equation may still 
incorrectly reproduce the two-particle distribution function. 

The inadequacy of Eq. (5.19) can also be understood intuitively as follows. Chaotic motions 
of any two particles experiencing the electrostatic pulses are “synchronized” since the pulses are 
independent of x and v. The Eangevin equation (5.19), on the other hand, desynchronizes particle 
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trajectories by involving additional Wiener processes, in spite of giving the correct one-particle 
Fokker-Planck equation. 

5.1.5 Example 2 

Next, consider a minority population of magnetized fast ions moving through a plane lower-hybrid 
wave that propagates perpendicular to the magnetic field. Assume the wave has a high harmonic 
number and a wavelength small compared to a typical ion gyroradius. Karney Kamey (1979) has 
shown that the dynamics of the perpendicular velocity of these ions are governed by a canonical 
time-dependent Hamiltonian system with Hamiltonian 

7/t = J — esin(\/^sin0 — z/f), (5.20) 

where / is the normalized magnetic moment, t the time normalized by the gyroperiod, 9 the gy- 
rophase, u the harmonic number, and e the normalized wave amplitude. Moreover, when e exceeds 
a threshold value, an ion’s motion becomes chaotic. This chaotic motion comes as the result of 
the effective randomization of the wave phase felt by an ion after a gyroperiod. Thus, above the 
threshold for chaos, we can model the wave phase as being randomized every gyroperiod by a 
random variable rj. That is, we can replace the exact chaotic ion motion with a stochastic ap¬ 
proximation; see Chirikov (1979) for Chirikov’s application of the same modeling approach to the 
standard map. This allows us to apply the formalism developed in this Chapter to find the physi¬ 
cal Langevin equation describing the stochastic particle trajectories at times much longer than the 
gyroperiod. 
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As in the previous example, the first step is to ealeulate si and S 2 . Set r = 27r and adopt the 
rough approximation 

°° 7 7 

^ exp(m6*) ^ ^ exp(mo6'), (5.21) 

p — n 0 

n=—oo 

where u = Uo + S, |5| < and = Jn{y/^) denotes the Bessel funetion of the first kind 
Abramowitz and Stegun (1964). This approximation amounts to selecting the most slowly varying 
term in the sum in Eq. (5.21). Then, upon directly evaluating the integrals in Eqs. (5.7a) and (5.7b), 
the resulting expressions for si and E[s 2 ] are 


Si = 27rsinc(7r(5) sin(no 6 * + rj) 


(5.22a) 




72 

^m+1 


72 

^m-1 


m=—00 
71 


m — p 


sinc( 27 r( 5 ) 


72 _ 72 

^rio+l ^Uo — l 


(5.22b) 


where 17 is a random variable uniformly distributed over the interval [ 0 , 27r] and sinc(x) = sin(x)/ x. 

Next, the space 77 can be constructed using the above expression for si. In this case, 77 is 
two-dimensional and has a basis {Xhi,Xh 2 }-, where 


(5.23a) 


= V27ismc{7i6)Jno{'/^) cos{no9) 
772(7,6') = \/ 27 rsinc( 7 r( 5 ) J„^(\/^) sin(no 6 '). 


(5.23b) 
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Finally, the coefficients for the Langevin equation, Eq. (5.12), can be derived using Eq. (5.13). 
The result is 


SI =e\/Tsinc(7r(5)noJno(v^) 

X (sin{no6)SW^ — cos{no6)SW‘^) 


(5.24a) 



(5.24b) 


The diffusion of the magnetic moment / predicted by Eq. (5.24) has already been studied by Karney 


Kamey (1979). However, Eq. (5.24) extends and compliments Karney’s results by predicting the 


appropriate diffusion in gyrophase, as well as the correct two-particle statistics. 

5.1.6 Concluding remarks 

We have shown how to derive the physical Eangevin equation for particle trajectories undergoing 
stochastic acceleration. This SDE correctly generates the correct one- and two-particle Eo kk er- 
Planck equations and inherits the Hamiltonian structure of the microscopic equations of motion. 
This inheritance is theoretically satisfying because it is a direct consequence of demanding con¬ 
sistency with the physical one- and two-particle Eokker-Planck equations. It also implies that 
symmetries of the macroscopic physical laws governing stochastic acceleration imply the presence 
of conservation laws. While this relationship is well known at the microscopic level, it is a pleasant 
surprise that it remains intact upon passing to dissipative macroscopic equations. 

A Hamiltonian Eangevin equation Eazaro-Camf and Ortega (2008) is a Stratonovich SDE of 
the form given in Eq. (5.12). If a loop of initial conditions for this SDE evolves under a given 
realization of the noise, then the action of that loop is constant in time. In addition, these equa- 
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lions arise from a stochastic action principle Lazaro-Camf and Ortega (2008) for which Noether’s 
theorem applies. Thus, by showing the physical Langevin equation is Hamiltonian, we have also 
identified potentially powerful tools for the analysis of stochastic acceleration. In particular, us¬ 
ing the methods of Bou-Rabee Bou-Rabee and Owhadi (2009), the stochastic action principle can 
be used to develop variational integrators for Eq. (5.12). Because these integrators are known to 
possess superior long-term statistical fidelity Bou-Rabee and Owhadi (2010), this approach may 
prove to be useful in Monte Carlo simulations of stochastic acceleration. 

5.1.7 Derivation of the physical two-particle Fokker-Planck equation 

This supplement to the article “The Hamiltonian mechanics of stochastic acceleration” consists of 
a derivation of the physical two-particle Fo kk er-Planck equation. In the main text, the two-particle 
Fokker-Planck equation is given in Eq. (6). The derivation will freely draw upon notation defined 
in the article. The essential idea behind this derivation is not novel; the same idea is presented in 
Bazant (2006) in the simpler context of a one-dimensional random walker. 

Fet Ft^s be the time advance map Abraham and Marsden (2008) associated with the dynamical 
vector field Xnt (Eq. (1) in our manuscript); Ft^s{z) gives the time t phase space location of a 
particle located at 2 ; G M at time s. Because this time advance map satisfies the identity Ft^r°Fr^s = 
Ft^s, where o denotes the composition of functions, we have F^^-^q = T)v.r,(A-i)T ° F(^n-i)t,{n- 2 )t ° 
... o Ft-q. This decomposition of the time advance map provides a stroboscopic description of 
particle dynamics; as the integer N increases, it tells us the phase space location of a particle at the 
times f = 0, f = r, f = 2r,.... 
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Because r -C n, the results in Lam (1998) may be used to write o 

exp(rX/i-(,), where St = esi^t + e^S 2 ,t + ... is given to O(e^) by 



(5.25) 


S 2 ,t = - / {exp{bXHo)*ht-b,ew{aXHo)*ht-a} dbda. 


Here exp(y) : M —>■ M is the time-one advance map of the dynamical system defined by the 
vector field Y; {•, •} denotes the Poisson bracket; and (exp(F)*h)(2;) = h(exp(—F)(2;)). This 
means that each r-second step in the discrete-time dynamics is a deterministic drift, exp{TXHo), 
followed by a small random kick, 6k = exp(Xs^j_^^j^). Moreover, the statistical assumptions on 
the perturbed force imply that the various 6k are independent identically distributed maps. There¬ 
fore, each discrete-time particle trajectory, znt = Fnt, oiz), is a temporally-homogeneous Markov 
process. 


In terms of Ft^g, the trajectory of a pair of particles at z = (^i, ^ 2 ) G M x M at time s is given 
by 


= {.Ft,s{,Zi),Ft^s{,Z2)). 


(5.26) 


This motivates introducing the two-particle time-advance map, Fg^t : M x M ^ M x M, which 
is defined by the formula 


Ft,s(-2) — {Ft^g{zi), Ft^g{z2)) ■ 


(5.27) 


The two-particle time-advance map inherits many of the qualitative features of the one-particle 
time-advance map. In particular, F*o = F^^^, which implies that F^t-^ = Fjvr,(iv-i)T o 
F(Ar_i)r,(Ar- 2 )r o ... o F.r,o- Morcovcr, F(^k+i)T,kT admits the decomposition F(^k+i)T,kT = ^k ° 
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exp(rXi^(,), where 


Sk{z) = {6k{zi),6k{z2)), 


(5.28) 


and bxp^tXho) is the time-r advanee map assoeiated with the veetor field on M x M given by 
the formula 


— Xhq{zi) © Xhq{z2) G © Tz^M. 


(5.29) 


The time-homogeneous Markov property implies that the linear operator defined on two- 
partiele observables Q : M x M —)■ M by the formula {P]\fQ){z) = E[Q(FAr.r,o(2;))] satisfies the 
semigroup property Pn+m = PnPm- Therefore, if we define the time evolution of a two-particle 
observable as Qnt = PnQ, then ^{Q(n+i)t — Qnt) = ^{Pi — ^)Qnt- For times sufficiently 
large compared with r, intuition suggests that the left-hand side of this identity approaches the 
partial time derivative dQt/dt and that ^{Pi — 1) may be approximated by a differential operator. 
This intuition can be made precise through the use of a Kramers-Moyal expansion Risken (1996); 
Bazant (2006), which we will describe now. Following Bazant (2006), we will obtain this limiting 
partial differential equation by scaling the time variable by the appropriate power of e and looking 
for a dominant balance of the equation 


1 


(Q{Ar+l)T — Qnt) 


1 


(-Pi — PjQnt 


(5.30) 


r 


r 


as e —)• 0. 
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First we eliminate the mean drift assoeiated with the baekground fields by working with the 
function 


Qt(z) = ^exp(-tXffJ*Qtj(z) = Qt(^exp(-tXffJ(z'^ 


(5.31) 


instead of Qt. In terms of Qt and t = Nt, Eq. (5.30) becomes 


Qt - exp^-rXnSQt + X] 


k=l 


9t'‘ 


T = 


(5.32) 


Qt - exp{-TXHj*Qt + E 


e[Qt, Si] + ^ 2 ] + [[Qt, Si], Si]/2) 


+ 0 (e^), 


where 


51 = exp^-tXnJ* (vTiSi,^ + tt^Si,^) (5.33) 

5 2 = exp^-tXnJ* (vriS2,r + T^2S2,r ); (5.34) 

the bracket [•, •] is the Poisson bracket on M x M defined by the formula 

[f,9]{zi,Z2) = {f{zi,-),g{zi,-)}{z 2 ) + {f{-,Z 2 ),g{-,Z 2 )}{zi); (5.35) 

and TTi, 7r2 : M x M ^ M are the projection maps onto the first and second factor respectively. The 
Taylor expansion in time is the key step here. It is justified by the fact that we will be considering 
late times when the evolution of Qt has had time to slow down as a result of diffusion. Notice that 
because E[si^t] = 0, E[si] = 0 as well. 
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Next we set t = (r/e^)A, where A is a renormalized dimensionless time. As e ^ 0, the 
dominant balance of Eq. (5.32) is given by 


dQx 

dX 


E 


[Qa, S 2 ] + [[Qa, Si], Si]/2 


(5.36) 


Or, in terms of t and Qt = ex.p{tXHo)*Qt, 

BO 

= [Qt,7^*lHo + e^E[7r]'S2,r]/T] + —E[[Qt,TTlSi^r],T^lSi^T]] 

.2 

+ [Qt,7T2Ho + e^E[n*S 2 ,r]/T] + —E[[Qi,7r;si,^],7r2Si,^]] 

+ ^E[[[Qt, 7r*Sl,r], TT^Si,^]] 

= L2Qt. (5.37) 


Thus, for late times Qt is given formally by 

Qt{z) = {exp{tL2)Q,6Q = {Q,gt,z) , (5.38) 

where {•, •) denotes the pairing of functions on M x M relative to the two-particle Liouville 
measure dzi dz 2 , 5z is a delta function concentrated at z, and gt^z is the distribution function of a 
pair of particles that begin dX z E M x M when t = 0. Because this identity holds for arbitrary 
functions Q, it implies that gt^z evolves according to 

gt,z = exp{tA2)6z, (5.39) 

where A 2 = L 2 is the adjoint of the operator L 2 . Differentiating this last identity in time 
and integrating against the initial two-particle distribution function finally leads to the two-particle 
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Fokker-Planck equation 


dgt 

dt 




(5.40) 


where A 2 is given by 

{A2gt){z) = A^^^gt + Af^gt + —E[a : (5.41) 

T 


The quantities Ai and a are defined in the main text. 

5.2 The Lorentz plasma 

5.2.1 Introduction 

The Lorentz plasma consists of a noninteracting gas of electrons moving through a neutralizing ran¬ 
dom collection of fixed, infinitely-massive, Debye-screened ions. The typical approach to studying 
the dynamics of the Lorentz plasma is to derive a Fo kk er-Planck equation governing the single¬ 
electron distribution function on time scales long compared with the plasma period. This approach 
eliminates the need to resolve the complicated structure of the ionic potential, and thereby greatly 
reduces the analytical and computational resources required to understand the plasma’s behavior. 

The structure of the Fo kk er-Planck equation, which takes the form of a Vlasov equation cor¬ 
rected by a collision operator, is ultimately determined by the microscopic Hamiltonian equations 
of motion for a single electron. Therefore qualitative features of the electronic equations of motion 
ought to have counterparts at the level of the collision operator. For instance, because an elec¬ 
tron that passes through the screened potential of an ion suffers no change in its kinetic energy, it 
would be surprising if the Fo kk er-Planck equation didn’t have a kinetic energy conservation law. 
Likewise, because there is no mechanism for electron absorption, the collision operator should be 
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consistent with the eonservation of eleetron number. The purpose of this Chapter is to study the im- 
plieations on the eollision operator of a more subtle qualitative feature of the electronic dynamical 
equations than either energy conservation of particle conservation. 

The qualitative feature we will be eoneemed with is the Hamiltonian nature of single-eleetron 
dynamies. Using the teehnique deseribed in Burby et al. (2013b), we will show that beeause the 
electronic equations of motion are Hamiltonian, there is a Fo kk er-Planck equation for the Lorentz 
plasma that is Hamiltonian in a stoehastic sense Lazaro-Cami and Ortega (2008). In partieular, this 
Fokker-Planck equation is the Kolmogorov forward equation associated with a stoehastic differen¬ 
tial equation that ean be derived from a stoehastic variational principle. We will then compare and 
eontrast this Hamiltonian Fokker-Planek equation with the elassieal result 

dtf + v-Vf = CL{f), (5.42) 

where Cl is the Lorentz eollision operator. We will find that the eollision operator, Chl, in the 
Hamiltonian Fokker-Planek equation is not identieal to but agrees with the latter asymptoti¬ 
cally in the limit Co, ei ^ 0, where e,, = Tac/r and ei = Vikt/L, t is the eoarse-graining time 
step, Tac ~ is the Lagrangian autoeorrelation time, and L is the length scale of the eleetron 
distribution funetion. We will also show that there is tension between the stochastie Hamiltonian 
property possessed by Chl and the energy-eonserving property of in the following sense. The 
operator Chl slowly produees energy for all non-zero Co, ei, while the limiting energy-eonserving 
operator Cl is provably not Hamiltonian in the sense of Lazaro-Cami and Ortega (2008). Finally, 
we will prove that any “reasonable” stoehastie Hamiltonian collision operator that is assoeiated 
with a path-wise energy-eonserving stoehastie differential equation must be signfieantly different 
from the Lorentz eollision operator Cl- 

Altogether, these results might give the impression that the stoehastie Hamiltonian formalism 
is not appropriate as an underlying mathematieal strueture for the piteh angle seattering proeess. 
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After all, sacrificing energy conservation is a hard pill to swallow. However, by carefully analyz¬ 
ing why the Hamiltonian Fokker-Planck equation for the Lorentz plasma slowly produces energy, 
we uncover a clue as to how stochastic Hamiltonian dynamics and energy conservation might be 
reconciled. We will argue that by relaxing the Markov assumption, and thereby allowing for some 
memory effects, the stochastic Hamiltonian approach may be able to accommodate energy con¬ 
servation. This possibility is interesting on mathematical, as well as physical grounds because 
there is presently no mathematical formalism for describing stochastic Hamiltonian dynamics with 
memory. 

5.2.2 Microscopic Hamiltonian description of the Lorentz plasma 

The set Q = will serve as the configuration space for the non-interacting electrons. The velocity 
phase space for a single electron is therefore M = TQ ^ Q x If 0 is the potential produced 
by the background ions, then the dynamics of each electron are governed by the Hamiltonian 

(5-43) 

2 rrie 

where qe and rrie are the electron charge and mass, respectively. The relationship between this 
Hamiltonian and the electron dynamical vector field X is given by Hamilton’s equations. 


ixUJo = dio, 


(5.44) 


where Uo = dx* A dvi. 

Structure of the ionic potential 

Let Ad and bo be the Debye length and the electron distance of closest approach, respectively. 
The plasma parameter A = Xn/bo- The structure of the electrostatic potential produced by each 
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background ion is determined by Ad and bo in the following manner. The charge density produced 
by an ion centered at Xj is given by p^. (x) = qi5{x — Xj). The electrostatic potential, , produced 
by such a Debye-shielded ion satisfies the differential equation 

= ^^Pxy (5.45) 

The only spherically-symmetric solution of this equation that decays as |a;| ^ cxd is given by 

0 ^.( 3 ;)= ^ exp (-lx - XjI/Ad) ■ (5.46) 

This “raw” potential, (f )^^, is not the mathematically-appropriate potential to subject to electrons in 
the Lorentz plasma. When an electron passes within a distance bo from the j’th ion, it experiences 
a large angle scattering event. Because such scattering events are exceedingly-rare, and because 
we would like to avoid infinities in our analysis, we will regularize the raw potential. We will 
accomplish this regularization by assuming that the potential produced by the j’th ion is given by 
^Xiix) = g{\x - Xj\) = - Xj\l\D), where 

g_{f) ifr<x 

^a(^) = S 4 if4<f<l (5.47) 

r A 

g^{r) iff>l, 

\ 

and g_,g_^_ are chosen so that (i) gj^{r) = 0 for f > 1 + <5 for some small <5 > 0 and (ii) the 
derivative of vanishes in a neighborhood of r = 0. The total electrostatic potential produced by 
N ions with centers Xj is then given by 

N 

<t^ = ^4>xy 
i=i 


(5.48) 
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Statistical properties of the ionic potential 

We will assume that the Xj are i.i.d. Q-valued random variables with PDF p. We will also work 
in the “thermodynamic limit.” The precise meaning of this statement is as follows. We allow the 
ionic PDF, p, to depend on the number of particles parametrically, i.e. p{x) = Pn{x). Then we 
make two assumptions: 

(i) iV > 1 

(ii) The limit 


lim Npi^ = Tii (5.49) 

N^OO 

exists (pointwise) and is equal to the constant n* = A/A||. 

The mean value of the ionic potential in the thermodynamic limit is given by 

(0) = lim E[(j)]{x) 

N^oo 

= \im N [ g{\x - x’\)pn{x’) dx^prime. 

N^oo J 

= „,/9(|x-x'|)dx' 

oo 

= 47 rnj J ( 5 . 50 ) 

0 

where dz denotes the standard volume form on It follows that the mean electrostatic force on 


an electron is zero. 
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The covariance function of the ionic potential in the thermodynamic limit is given by 

C^<p{x,y) = lim E[(0 - E[0])(a;)(0 - E[0])(|/)] 

A'^oo 

= nij gi\x - x'\)g{\y - x'\)dx'. ( 5 . 51 ) 

where we have used the fact that the Xi are independent and E[(;/)a.J oc N~^. Notice that 

C4>4>i^,y) = C{\x-y\), ( 5 . 52 ) 

where 

C{d) = ^^C^{d/\D), ( 5 . 53 ) 

Ztt ad 

and 

CAd) = 

oo oo _ 

I I ^ a (\/^^ + + ( 5 - 54 ) 

—OO 0 

is a dimensionless covariance function. The most important properties of and C are the 
following. 

(PI) is manifestly positive semi-definite, i.e. for any finite collection of points Xj G Q 
and corresponding real numbers a^, 

EE 

i 3 


( 5 . 55 ) 
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Thus, is an isotropic covariance function, and C is a radial basis function. 

(P2) C{d) = 0 when d > 2Ad+- This follows from the fact that g{i-) is compactly sup¬ 
ported in the interval [0, Ad+)- 

(P3) C{d) = C{-d) 

Using (P3), we can also see that all odd powers of d vanish in C’s Maclaurin series, i.e. 

C{d) = C'(O) + ^^"(O) d^ + ^C"'"(0) d^ + 0(d®), (5.56) 

as d —)■ 0. Thus, 

-C"'(d) = 0(^2), (5.57) 

and 

C'{d) = C”{0)d + O{d^), (5.58) 


as d —)■ 0. 

The covariance tensor of the electrostatic yip/J produced by the ions is given by 


^v^vAx,y) = lim E [V(0 - E[0])(a;)V(0 - E[0])(|/)] 

N^oo 

= lim E [V0(a;) V0(i/)] 




(5.59) 
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where 


C(d) 




(5.60) 


and 


U{d) 


C’A{\d\) / d d\ 

\d\ V |c^l|c^|/ 


cli\d\) 



The most important properties of Cv</-V(/. and <t are the following. 


(5.61) 


(FI) Cv(;iv</- is positive semi-definite, i.e. given any finite eolleetion of points Xj and eorre- 
sponding veetors Vj, 


EE Vj ■ Xfc) • Lfc > 0. (5.62) 

j k 


(F2) C(A) = 0 when |A| > 2 Ad+- 

(F3) C(A) = C(-A). 

Using (F2) and integration by parts, we ean simplify the following type of definite integrals 
involving C. Let L > 2A£)+ and ehoose a unit veetor e. Set 

L 

4(e) = j |ArC(Ae) d\. 


-L 


(5.63) 
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For integer n > 0, we have 


hie) = -2 



n—1 


C'{X)dX (id — (n + l)ee). 


(5.64) 


5.2.3 The argument for employing a Fokker-Planck equation to model the 
Lorentz plasma 

The purpose of this section is to give a qualitative motivational picture of the ideas that go into 
deriving a Hamiltonian Fo kk er-Planck equation for the Lorentz plasma, as well as to describe pre¬ 
cisely what a Hamiltonian Fo kk er-Planck equation is. Asa point of departure, we remind the reader 
of the justification and derivation of a Fo kk er-Planck model in terms of so-called jump moments. 
We then reformulate the same discussion in terms of the symplectomorphism group, Difft^^(T(5). 
This reformulation offers a particularly suggestive description of how the Hamiltonian nature of 
the electronic equations of motion influences the structure of the Fo kk er-Planck equation. 


In terms of jump moments 

The problem of finding a formula for the trajectory of an electron in the Lorentz plasma is ex¬ 
tremely complicated. There is not just a single ion; if there were, we would only have to solve the 
two-body problem. There are not just two ions either; if this were true, we would be faced with 
(an analogue of) the soluble Euler three-body problem. It is better to assume there are ~ 10^^ ions, 
and therefore electron dynamics are surely chaotic. As such, when studying the dynamics of the 
Lorentz plasma, we must be satisfied with less detailed information than exact electron orbits. 

One way to give a less detailed (and therefore simpler) description of the Lorentz plasma that 
still retains a great deal of dynamical information is to find the evolution equation for the single¬ 
electron PDF. This is a much more manageable task than finding the precise electron trajectories 
for the following reason. Whereas the large number of ions mangles the electron trajectories, it ac- 
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tually simplifies certain statistical properties of electron dynamics. For instance, the (Lagrangian) 
autocorrelation time of the force on an electron is set by the plasma period t^c ^ Thus, if 
we chop the time-axis into intervals of length r 3> Tac and restrict attention to an electron’s phase 
space location at the ends of these intervals, we obtain a discrete-time Markov process in phase 
space, Zn, where n is the discrete time index. If we could estimate the transition probabilities 
of this Markov process, we would be able to derive an equation describing the evolution of the 
electron PDF on time scales much longer than r. This derivation could proceed by analogy with 
the derivation of the late-time evolution equation for the PDF of a random walker given in Bazant 
(2006). Indeed, our electron Markov process is nothing more than a random walk in phase space 
where the PDF of the walker’s step depends on where the walker is standing. 

Provided that we make the further restriction r -C Tbounce, where Xbounce is the characteristic 
time for an electron trajectory to deviate significantly from free streaming, the transition probabil¬ 
ities can be calculated using ordinary perturbation theory * . It is straightforward to show that the 
resulting late-time evolution equation for the single-electron PDF takes the form 

dtf + diw{fu) = div(T> ■ d/), (5.65) 

where div denotes the divergence relative to the Liouville volume form dx dv, m is a vector field 
on phase space, and D is a rank-2 tensor on phase space with components . The drift vector u 
is given by 


U = Uo+ {A2 z)/t 


^{di\{Aiz)Aiz)/T, 


(5.66) 


*Our ability to require that r satisfies both r <C Tbounce and r ^ Tac follows from the fact that small-angle scattering 
events dominate over large-angle scattering events. 
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where Ai^, A 22 ; are veetor fields on phase spaee and Uo = v ■ dx is, the free-streaming veetor field. 
The diffusion tensor is given by 


D = ^(Aiz (g) Ai2:)/r. (5.67) 

We will refer to Aiz and A 22 : as the jump vectors. They are defined so that a particle starting at 
Zo = (x, v) when t = 0 ends at 


Zr ^ +/S.2z){x + vT,v) (5.68) 

when t = T, with second order accuracy^. The components of D are known as the jump moments. 
Equation (5.65) is known as the Fokker-Planck equation. 

In terms of the symplectomorphism group 

The previous argument justifying the use of a Fokker-Planck equation for the late-time single¬ 
electron PDF is appealing because it draws upon only elementary facts about Markov processes. 
However, it has the disadvantage of obscuring a striking geometric picture that underlies the whole 
discussion. In order to illuminate the geometric picture, we will now give a second justification for 
using a Fo kk er-Planck equation. 

Fet be the t = r time-advance map associated with the full single-electron equations of 
motion. Given an initial condition 2 ; G TQ, the approximately-Markov process that gives an 
electron’s phase space location at the ends of r-second time intervals is given by 

Zn = T^{z), (5.69) 

^The amplitude of the fluctuating electric field can be regarded as the expansion parameter, which is essentially 
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where n G Z is the diserete time variable and denotes the n-fold eomposition of T with itself 
(e.g. when n = 2, Note that the mapping is random beeause it depends on 

the eonfiguration of the ions. Because Zn is a Markov process, the operator Un given by 

{UnQ){z) = (Q(^”(z))), (5.70) 

where Q : TQ —)■ M is an arbitrary observable, must satisfy the semi-group property 

Uni+n2 (5.71) 

However, by noting Ui = where J^* denotes the pullback operator along 7^, we also have 

U2 = {rr) ^ = UiUi. ( 5 . 72 ) 

Note that T and T are not independent random mappings - they are identical! This contradiction 
tells us that the process Eq. (5.69) is not precisely Markov. 

In spite of this contradiction, we know that Zn is approximately Markov because r 3> Tac. 
Therefore, for the sake of modeling it is sensible to replace Eq. (5.69) with 

Zr, = Tn^y^n-\°---^^\{.z), (5.73) 

where o denotes functional composition and the Ti are i.i.d. random mappings each with the same 
PDE as T. The effect of this replacement is that the locations of all of the ions are scrambled 
after each time step. While this scrambling effect is, strictly speaking, unphysical, it ought to be 
statistically harmless; electrons forget about the orientations of the ions after Tac seconds anyway. 
It is easy to check that this redefined Zn is rigorously Markov. 


^This formula shows that Ui is a mean propagator. 
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Equation (5.73) tells us that one way to determine evolution of the single-eleetron PDF is to 
first determine the statistical behavior of the mapping-valued stochastic process 

9n ^ n ^ ^ n—1 O • • • O ^ (5.74) 

Indeed, the expected value of Q{zn), for any observable Q, is given by 

(Q(z„)) = (s;Q)(j) = Wn)Q)(~-), (5.75) 

which shows that the single-electron PDF is completely determined if we know But why 
would we want to do this? Whereas the process Zn evolves in a finite-dimensional space, the 
process gn evolves in an infinite-dimensional space of mappings, which suggests that Qn is a much 
more complicated object than Zn- The answer is that gn is simpler than Zn when viewed in the right 
way. 

First notice that gn is a diffeomorphism for all n, i.e. gn is smooth and invertible with a smooth 
inverse. This follows from the fact that J^i is a diffeomorphism for each i (being a time-advance 
map for an ODE on phase space) and gn is a composition of the J^’s. This means that the mapping¬ 
valued process gn takes place in a very special space of mappings known as the phase space dif¬ 
feomorphism group Diff(TQ). Diff(TQ), which is the set of all diffeomorphisms of the velocity 
phase space TQ, is a group under functional composition. In a sense that we will not discuss here, 
Diff(r(5) is also a smooth (infinite-dimensional) manifold. Thus, the process gn evolves in a space 
with a very rich structure. 

Next notice that the increments 5gn2,ni = ° ^ have the following simple 

statistical properties: 


(RWl) If Hi < n 2 < ris < n^, 5gn2,ni and Sgn^^ns are statistically independent. 
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(RW2) The PDF'^ of Sgn 2 ,ni depends on ni,n 2 only through the differenee n 2 — ni. 

If we were to replace Diff(TQ) with the additive group M, these properties would imply that 
gn is an ordinary random walk^. More generally, if Diff(TQ) is replaced by any group G, a 
discrete-time process on G that satisfies properties (RWl) and (RW2) is known as a random walk 
on G. Thus, the process gn is a Diff(rQ)-valued random walk. 

Once nice feature of thinking about the dynamics of the Lorentz plasma as a random walk 
on the diffeomorphism group is that we have a good intuitive understanding of the long-time 
behavior of random walks. In particular, we know that, under an appropriate scaling limit, an 
ordinary random walk is well-approximated by a Brownian motion. Therefore we can reasonably 
expect that the long-time behavior of the random walk gn is described by a Brownian motion 
on Diff(TQ)Baxendale (1984), i.e. a continuous-time process gt G Diff(T(5) that satisfies the 
properties 

(BMl) If ti <t 2 <H < ti, 5gt2,ti and 5gt^,t3, are statistically-independent. 

(BM2) The PDF of 5gt2,ti only depends on ti,t 2 through the difference ^2 — h- 

(BM3) The sample paths of gt are almost surely continuous functions of t. 

Of course, properties (BMl) and (BM2) are the obvious analogues of the properties (RWl) 
and (RW2) that we already know gn satisfies. Property (BM3) is motivated by the dominance 

"^Here PDF stands for probability distribution/wncf/onaZ. 

^The PDF of a step taken by such a walker would be arbitrary. 
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of small-angle scattering events over large-angle scattering events, in particular the fact that the 
particle jumps can be calculated using perturbation theory. 

In Baxendale (1984), Baxendale shows that Brownian motions on Diff(TQ) always arise as 
the stochastic time-advance maps of stochastic differential equations. If gt is such a stochastic 
time-advance map, then a basic fact from the theory of stochastic differential equations states 
that the PDF of the random variable gt{z), where 2 : G TQ is a fixed initial condition, satisfies a 
Fokker-Planck equation of the form given in Eq. (5.65). We have therefore arrived at an alternative 
justification for the use of a Fokker-Planck equation to model the Lorentz plasma. 

Another nice feature of thinking in terms of gn instead of Zn is that the Hamiltonian nature of 
the electron dynamical equations manifests itself in a very simple way at the level of Diff(TQ). 
Indeed, because the electron dynamical vector field X given by Eq. (5.44) is Hamiltonian, the 
t = T time-advance map X must preserve the symplectic form, = ce,,, where ujo = dx'' A dvi. 
Eikewise, because the Xi have the same PDE as X, we also have X*u:o = uio- Therefore, by the 
identity {Xi o Xj)* = the process gn satisfies 


(5.76) 

for all n. In other words, gn is not free to wander everywhere in Diff(T(5), but only along the 
constraint set defined by g*ijJo = ^o- Actually, the latter constraint set is a subgroup of Diff(TQ) 
known as the symplectomorphism group, Diff^^(TQ) c Diff(TQ). Because gn does not leave the 
symplectomorphism group, the limiting Brownian motion gt also must satisfy the same constraint. 
It can be shown that the latter requirement constrains the Eokker-Planck equation (5.65) to satisfy 

u = (5.77) 

D = '^Xh^ ®Xh^, 

k=l 


(5.78) 
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where Xf denotes the Hamiltonian veetor field with Hamiltonian /, the hk are arbitrary funetions 
on phase spaee, and the sum over k is possibly infinite. Conversely, given a Fokker-Planek equation 
whose drift veetor u and diffusion tensor D take the above form, it is always possible to find a 
Diff^^(TQ)-valued Brownian motion that generates it. 

When the drift veetor and diffusion tensor of a Fokker-Planek equation are in the form pre- 
seribed by Eqs.(5.77) and (5.78), we will say that the Fokker-Planek equation is Hamiltonian. 
Thus, an important consequenee of the faet that the single-eleetron equations of motion are Hamil¬ 
tonian in nature is that the Fo kk er-Planek equation for the single-eleetron PDF ought to be Hamilto¬ 
nian. By making this observation, we ean see that in passing from the mieroseopic single-electron 
equations of motion to the macroscopic Fokker-Planek equation, we move from the world of 
Hamiltonian mechanics into the world of stochastic Hamiltonian mechanics. Where Hamiltonian 
mechanics is concerned with smooth one-parameter subgroups of Difftj^(T(5), stochastic Hamil¬ 
tonian mechanics is concerned with Brownian motion on Diff^^(T(5). The stochastic Hamiltonian 
nature of the Fo kk er-Planek equation is the moral counterpart to the Hamiltonian nature of the 
microscopic equations of motion alluded to in the introduction. 

5.2.4 Hamiltonian Fokker-Planek equation for the Lorentz plasma 

We will now apply the technique described in Burby et al. (2013b), which we will refer to here¬ 
after as BZQ, to derive a Hamiltonian Fokker-Planek equation for the Lorentz plasma. In BZQ’s 
notation, we have 



(5.79) 


(5.80) 
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where cj) is the ionie potential deseribed earlier. The first- and seeond-order jump veetors are given 
by Aiz = Xsi and A 2 Z = Xg^, where 

51 = [ Ff^hdt (5.81) 

Jo 

52 = ^ [ [ {Ft^^h,Ft^^h} dt2dti. (5.82) 

^ Jo Jo 

The unperturbed flow map, Ft, is given by 

Ft{x,v) = {x + vt,v). (5.83) 

Using these formulae, the drift veetor uhl and the diffusion tensor Dhl can be eomputed in terms 
of the eovarianee tensor € given in Eq. (5.60). 

First we eomputew^L = Mq + 36(^2)/r, whieh amounts to eomputing ( 52 )- It is straightforward 
to verify that the mean of the Poisson braeket appearing in the definition of S 2 reduces to 

{{Ft,,h,Ft,M) = -(—) [ti-t2H(t{[F-t2]v)), (5.84) 

yTfle J 

where tr denotes the trace. Upon substituting this expression into the double integral, changing 
integration variables, and applying Fubini’s theorem, we then arrive at the following expression 
for (S 2 ): 

(s 2 ) = -^f—) [ t{T-t)ti{(E{vt))dt. (5.85) 

2\meJ Jo 

Note that because ( 52 ) only depends on (x, v) through v, the drift vector uhl = Mo + {^ 2 z)/t = 
Uo + 3f(s2 }/t only has an x-component. 

Next we derive an expression for the diffusion tensor Dhl- For this purpose, we introduce a 
useful notation for contravariant second rank tensors on phase space. If T is a dyad like C, then 
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is the second rank contravariant tensor on phase space given by 


dhi ■ Txv ■ dh2 = Vhi ■ T ■ 


(5.86) 


where hi, h 2 are arbitrary functions on phase space. The tensors T^x-, T^x-, Tw are similarly defined. 
In terms of this notation, 0 X^^) is given by 


(X,,0X,,)= ^ 


m. 



^vv{[t 2 - ti]v) dti dt 2 


0 0 



+ / t 2 <txx{[t 2 - tl]v) dti dt 2 + ti(txx{[t 2 - tl]v) dti dt 2 



0 0 


0 0 


r r 

+yy tit2(txx{[t2-ti]v)dtidt^. 
0 0 


(5.87) 


After simplifying the double integrals and dividing by 2r, the diffusion tensor D is then given by 


Dhl = 


1 / 7e 


T xrrif 


(r — t)<tvx{vt) dt 


T 


(r - t){(txv{vt) + (txx{vt)) dt 


2 PT 


(r - 3i/2 + iy 2 T^)(txx{vt) dt 


(5.88) 


As explained in Burby et al. (20 1 3b), because Dhl = ® X^^)/(2r), there exists an expansion 

of D in the form 


Dhl — ^ Xhk ® (5.89) 

k=l 

where the form an orthonormal basis for the reproducing kernel Hilbert space associated with 
the phase space covariance kernel a{zi, Z 2 ) = (X5j(2;i) 0 Xs2{z2)) / {2t). We will not find the hk 
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here. It is worth pointing out, however, that finding the hk is equivalent to the well-known problem 
of finding the so-called Karhunen-Loeve Lord et al. (2014) expansion of a Gaussian random field 
with covariance a. 

Because uhl and Dh^ just calculated can be written in the form given in Eq. (5.77), the Fokker- 
Planck equation 


dtf + div{fuHL) = div{DHL ■ df) 


(5.90) 


is Hamiltonian. In the following section, we will compare and contrast this Fokker-Planck equation 
with the classical result given in Eq. (5.42). In particular, we will compare the Forentz collision 
operator 


CL{f)=diw{DL-df), 


(5.91) 


where 


Dl = z/(u) U^^(u), 


(5.92) 


3 

U(u) = |up(id — vv), and iy{v) = with the Hamiltonian collision operator 


CnLif) = div{DHL ■ df), 


(5.93) 


where Dhl is given by Eq. (5.88). 
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5.2.5 Comparison of the Hamiltonian collision operator and the Lorentz op¬ 
erator 

Asymptotic equivalence 

In order to demonstrate the asymptotie equivalenee of Eqs. (5.42) and (5.90), we introduee the 
dimensionless variables x, n, t. These normalized position, veloeity, and time variables are related 
to their unnormalized eounterparts by 


X = Lx 

(5.94) 

V = VthV 

(5.95) 

t = Tt, 

(5.96) 


where L is the length seale of the eleetron distribution funetion, Vth = Ais the thermal veloeity, 
and T is the temporal seale of the eleetron distribution funetion. We will set T = A/ujp, whieh 
is eonsistent with measuring time in units of the eleetron-ion eollision period. The Hamiltonian 
Fokker-Planek equation expressed in these dimensionless variables is 

dtf + div(/M) = div(T> ■ df), (5.97) 

where u = Tuhl and D = TDhl- We will now present expressions for u and D that exhibit their 
dependenee on the small parameters €« = l/{TUJp), ei = Vthr/L, and 1/A. 

The normalized drift veetor is given by 


1 

3 ^ 


u = Aeoeiv ■ - 


1V U A ’ ) 


(5.98) 
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where the dimensionless function x is given by 


_ _ _ _ _ 
x{v) = / eoA(l-eoA)trCA(wA)(iA. 
Jo 


(5.99) 


Note that x depends on Co and A, but not ei. 
The normalized diffusion tensor is given by 


JJ 10^2 2 T "^xv) T ^\C,xx ) ) 


(5.100) 


where the dimensionless dyads 77 , are given by 


7 ) 

c 


-l/co 


(1 — CoA) Ca(lA) d\ 



eoA /2 + /o)^\{v\) d\. 


(5.101) 

(5.102) 


Note that rj, C, like x, depend on So and A, but not ei. 

These expressions for u and D can now be used to study the asymptotic behavior of the Hamil¬ 
tonian Fokker-Planck equation as € 0 , ei and 1/A tend to zero. For concreteness, we will study this 
limit under the assumption 


-/Va 

(5.103) 

-/Va. 

(5.104) 


Effectively, this assumption chooses a specific path to zero through (cq, ei, 1/A)-space along which 
our asymptotic limit is taken. 
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First we consider the limiting behavior of the quantities x, Vi C- For this purpose, it is enough 
to consider the limiting behavior of the integrals 


/“V'o _ _ _ _ 

ln= / (e,ArCA(FA)dA 
Jo 

,n /‘l«l/eo 


Jo 


€a{vs) ds, 


(5.105) 


for non-negative integer n. It is not difficult to show that 


Irt, OC. 


In A if n = 0 


e” if n > 0. 


Therefore the asymptotic limits of x, i?, C given by 


(5.106) 


A ^ 0 

1 ^ 

77 ^ J = — / ^Af-sn) ds 

Jo 


(5.107) 

(5.108) 

(5.109) 


The limiting drift vector and diffusion tensor are now simple to obtain. For u we have 


u ^ V ■ dx- 


(5.110) 


Similarly, D is given by 


(5.111) 

The dyad I can be simplified further by using the fact that when A is large, C^id) ~ 27re“l‘^l 
for \d\ > A~^. In fact, if the Debye screened potential was not regularized, this would not be an 
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approximation. Thus, 


I 



27 rln A 


(id — vv). 


(5.112) 


V 


After restoring units, we can therefore write the limiting drift vector and diffusion tensor as 


u = V ■ dx 


(5.113) 



(5.114) 


where the dyad U = |Lp(id — vv). Comparing these expressions with Eq. (5.91) reveals that we 
have indeed recovered the classical result for the Lorentz plasma Fokker-Planck equation. 

This result assumes the scalings given in Eqs. (5.103) and (5.104), but many other choices seem 
allowable. Therefore a natural question is whether or not the form of the limiting Fokker-Planck 
equation depends on how we send the small parameters to zero. It turns out that there are only 
two possible limiting equations, and these differ only in the limiting drift vector; either the free 
streaming term survives the limit or it doesn’t. We chose our scaling so that the free streaming 
term survives because this seems to be the most interesting possible scenario. 

The Hamiltonian collision operator slowly produces energy 

While the Fokker-Planck equation (5.90) is manifestly Hamiltonian, it does not conserve kinetic 
energy. This can be seen by direct calculation. The total electron kinetic energy is given by 



(5.115) 
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where Vt = dxdv is the Liouville volume form. The rate of ehange of the plasma kinetie energy is 
therefore 



HodtfQ 

Hodiw{-fuHL + DHL-df)Q 


j dHo-DuL-dfn 


= / div{DHL-dH,)fn, 


(5.116) 


where we have used the faet that ( 52 ) only depends on v and Dhl is symmetrie. This expression 
shows that the only way kinetie energy will be eonserved regardless of initial conditions is if 
dii\{DHL ■ dHo) = 0. However, 


div(DHL ■dHo)=[ — \ V, 

rUp 


(1 — t/r)n • ^{vt) dt 


TUf 


Qe 


1 


1 

_ 0 1 


r 


d\v\ 



(1 — t/ t)vC''{\ v\t) dt 


'ti|[c( 0 ) -C(|K|r); 


C(0) 


nip J \vrT 


(5.117) 


where the last line is valid when |n|r > 2AD+. It follows that ^ 7 ^ 0 and that the eharaeteristie 
time for energy ehange is 


Te 


|g| 

\dE/dt\ 


tA. 


(5.118) 


We will have more to say about why Chl does not eonserve energy in the final seetion. 
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The Lorentz operator is not Hamiltonian 

While Chl does not eonserve energy, the limiting eollision operator Cl eertainly does. Therefore 
it is tempting to hope that the limiting proeedure that transforms Chl into Cl preserves the Hamil¬ 
tonian nature of Chl- In this seetion we will prove definitively that this is not the ease. Speeifieally 
we will show that there is no sequenee of funetions hk sueh that Dl = ® Xh^. 

The level of mathematieal diseourse inereases substantially in this seetion for two reasons. 
First, a basie knowledge of veetor bundles and operations on veetor bundles is assumed. A partie- 
ularly readable aeeount of this material is given in Bott and Tu (1982), starting on p. 53. Seeond, 
we assume a working knowledge of the theory of reprodueing kernels. Reprodueing kernels are 
really nothing more than two-point eovarianee funetions, but mueh ean be said about them with¬ 
out referring to ideas from probability theory, and this (somewhat) justifies introdueing a seeond 
name for them. The standard referenee for learning about the basies of reprodueing kernel theory 
is Aronszajn (1950). 

Suppose that M is a smooth manifold and k:MxM— )-Misa smooth reprodueing kernel. 
Assoeiated to k is a smooth seetion, a^, of the symmetrie tensor produet T*MqT*M. The seetion 
Ok is defined as follows. Let vi,V 2 G T^M be tangent veetors at x e M. Choose smooth eurves 
ci,C 2 : / M, where / is an open interval of M eontaining 0, sueh that ci(0) = C2(0) = x, 
c)^(O) = vi and = V 2 . We set 


OtKivi,V2) 


d d 

dei Qde2 


K(ci(ei), 02(62)). 
0 


(5.119) 


A seetion a of T* M Q T* M has a reproducing kernel primitive if there is some reprodueing 
kernel n sueh that a = The purpose of this seetion is to prove a theorem that eharaeterizes 
the set of a’s with reprodueing kernel primitives. We will only eonsider seetions a with locally 
constant rank. The rank of a seetion a at x G M is defined as the eodimension of the kernel of ax. 
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More precisely, if we let T^xM be the collection of vectors Vx such that 

'iwx e TxM, a{yx,Wx) = 0, (5.120) 

then the rank of ax is rank(Q; 2 ;) = dim{TxM) — dim{TaxM). A section a has locally constant 
rank if for each x E M there is some open neighborhood of x on which the function u E M 
rank(Q;„) G Z is constant. 

First we will prove that if has locally-constant rank, then Lxa = 0 for any vector field 
X that takes values in a^’s characteristic distribution. Given an a with locally-constant rank, its 
characteristic distribution is the subbundle of TM whose fiber at a; G M is given by TaxM. a’s 
characteristic distribution will be denoted Tq,M. A general a with locally-constant rank need not 
satisfy Lxa = 0. Thus, the following theorem provides a non-trivial necessary condition for a 
section a with locally-constant rank to have a reproducing kernel primitive. 

Theorem 1. Let nbe a smooth reproducing kernel on M. a^ is positive semi-definite as a bilinear 
form. Moreover, ifa^^ has locally constant rank, then Lxa,., = Ofor each vector field X that takes 
values in T^^M. In particular T^M is integrable in the sense of Frobenius. 

Proof. Let {fj} be an orthonormal basis for the reproducing kernel Hilbert space associated with 
K. The kernel k can be expressed in terms of these basis elements as k{x, y) = (j)j{x)(l)j{y), 

which implies that a„ = ’^jidfjfi. This immediately implies that is positive semi-definite as 
a bilinear form. 

Select anxo ^ M and restrict attention to an open neighborhood U of Xo where u i-A rank(Q;K«) 
takes the constant value r. Suppose n is the largest integer with the property that there exists n 
distinct basis elements f^ = ..., that satisfy 


{d(l)jf)xa A • • ■ A {d(l)j„)xa 7^ 0. 


(5.121) 
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Because the space of n-forms with n > m = dim(M) vanishes, n must satisfy n < m. Therefore 
we must be able to find m — n additional functions so that the /* comprise a coor¬ 

dinate chart on U (it may be the case that U must be shrunk). In this special coordinate system on 
U, Eq. (5.121) implies that each must be independent of the coordinates ..., /"*. Thus 
each of the m — n vector fields dm+i = ffm+i , ■ ■ ■ idn = -§fK take values in the characteristic 
distribution of a^- It follows that the dimension of the characteristic distribution, m — r, must be 
greater than or equal to m — n, i.e. n > r. In fact, n cannot be greater than r. To see this, note that 
Ok = tto + ^ 01 , where 

n 

«o = 5^(d/T (5.122) 

i=l 

5a = (5.123) 

It is straightforward to show that rank(ao) = n. Moreover, because ao and ba are each positive 
semi-definite bilinear forms (being sums of squared 1-forms), r = rank(a,,) > rank(ao) = n. 
Therefore we must have n = r. This shows that dr+i-, ■ ■ ■ ,dm in fact span a^’s characteristic 
distribution. Moreover, being pushforwards of some of the standard basis vectors in these 
m — r vector fields commute. We have therefore succeeded in proving that the commutator of 
any pair of vector fields that take values in Ta^M also takes values in this is precisely 

integrability in the sense of Frobenius. Actually, we have nearly proved more than this. If X is any 
vector field that takes values in T^^M, then it must be a C°°( (7) -linear combination of the vectors 
dr+i,... ,dm- Therefore, 


LxcTk = 2YXLxd(j)j){d4>j) 

3 

= 2Y{dLx<Pj){d(j)j) = 0, 

3 


(5.124) 
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where, on the last line, we have used the faet that eaeh (j)j is independent of , /™. In other 

words, for eaeh X that takes values in T^^M, Lxcn^ = 0. It is not hard to show that this last 
property aetually imples T^M is integrable in the sense of Frobenius. 


□ 


Next we will show that the necessary conditions for a constant-rank a to have a reproducing 
kernel primitive given in Theorem 1 are, in a particular sense, locally sufficient. 

Theorem!. Suppose a G T{T*MqT*M) ispositive-semidefinite, has locally-constantrank, and 
that Lx<y = Ofor each vector field X that takes values in TaM. Then for each Xo G M, there is an 
open set U containing Xo such that a\U = a^jj, where nu '■ U x U ^ is a smooth reproducing 
kernel on the open submanifold U. 

Proof. Let r be the rank of a near Xo- By the Frobenius theorem we can choose coordinates /* on 
a neighborhood U of Xo with the following property. The vector fields ..., ^ do not lie in 
TaM, while the vector fields dr+i = gyr+r, ■ ■ ■ idm = do. In this coordinate system, a must 
be of the form 


r 



(5.125) 


where Aij is an r x r symmetric, positive definite matrix of functions on U (note that the upper limit 
of the double sum in Eq. (5.125) is r < m). By assumption, LQ^a = Q for each I G {r +1,..., m}, 
which implies 



(5.126) 
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Because the {dp){df^) with i < j are linearly independent and the matrix Aij is symmetric, the 
previous equation shows that 


dAij 

Of 


0 , 


(5.127) 


for each r} and I G {r+1, ..., m}. i.e. the Aij only depend on the first r coordinates 

in this coordinate system, Aij = Aij{f ^,..., 

Let V C be the image of the submersion it : u ^ U G M’’. Without 

loss of generality, we can assume that V is open and connected. Because the Aij only depend on 
the first r coordinates on U, they define a positive-definite bilinear form g on V given by 


g{x^, ..., x^) = Aij{x^, ..., x^)dx'' dx^. (5.128) 

In other words, {V, g) is a Riemannian manifold with metric tensor g. By Nash’s embedding 
theorem, there is therefore an isometric embedding I : {V,g) (K^,go) where M is some 

integer and go = is the standard metric tensor on In particular, 

M 

g = rgo = J2{ds,f, (5.129) 

i=l 

where Si = I*x\ 

Now, it is simple to verify that a = 7r*g, and therefore 

M M 

a = 71* ^(ds^^ = (5.130) 

i=\ i=l 

where 0* = 7i*I*x\ We have therefore proved that a\U = where k,{x, y) = 
is a reproducing kernel. 


□ 
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If we assume that the leaf space of the foliation tangent to a’s characteristic distribution is a 
manifold, the previous theorem can also be globalized in the following manner. 

Theorem 3. Suppose a G r(T*M 0 T*M) is positive semi-definite, has locally-constant rank, 
and that Lxa = 0 for each vector field X that takes values in T^M. Also assume that the leaf 
space of the foliation tangent to T^M is a smooth manifold with the quotient topology. Then there 
is a reproducing kernel k : M x M —)■ R such that a = a^. 

Proof Let Mq be the leaf space of the foliation tangent to a’s characteristic distribution. Let 
TT : M —)■ Ma be the map that sends a point in M to its corresponding leaf. Because Lxoi = 0 for 
each X taking values in T^M, there is a unique g G V{T{Ma) O T{Ma)) such that a = TT*g. 

To see that g is uniquely determined by the formula a = Ti*g, consider the following. Let 
wi,W2 G T^Ma be a pair of vectors tangent to the leaf space at a; G Mq. Let wi,W2 G T^M be any 
pair of vectors tangent to M at a; G M that satisfy TTT{wi) = Wi. Suppose w'2 G T^'M are also 
lifts of the vectors Wi,W 2 . Then there must be a diffeomorphism ^ : M ^ M that preserves the 
leaves of the foliation tangent to TaM (i.e. tt o $ = tt) such that x' = <f>(a:). This diffeomorphism 
allows us to compare the vectors Wi, W2 and wf wf In particular, we can consider the differences 
(5i = w[ — T^{wi) and 82 = w^— T^{wi). We have 

T7r(<5,) = T7t[w' - T^Wi)] 

= Wi- T{tt o 


= Wi — Wi = 0. 


(5.131) 
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This identity allows us to compute the difference between a{w[, W 2 ) and a{wi, W 2 ). Indeed, 

a{w[,W2) = a{ 5 i + T^{wi), 82 + T^{w2)) 

= a{T^{wi),T^{w2)) 

= ^*a{wi,W2) 

= a{wi,W2), (5.132) 

where the last line follows from Lxa = 0. It follows that a(wi, W 2 ) depends only on a, wi, and 
W 2 . Thus, g{wi, W 2 ) = a{wi, W 2 ) is well-defined. 

The bilinear form g is positive definite for if w G TM„ is tangent to the leaf space, g{w,w) = 

a{w,w), and the right-hand-side is zero only when w is tangent to T^M, i.e. when w = 0. 

The pair {Ma,g) is therefore a Riemannian manifold. Nash’s theorem then implies that there 
exists an isometric embedding / : (Mq,, g) go) where go is the standard euclidean metric 

9o = . We have therefore proved that 

No 

a = Ti*g = 7T*rgo = (5.133) 

i=l 

where (pj = tt* I* xK Equivalently, a = where the kernel k 

We will now use Theorem 1 to prove that the Lorentz diffusion tensor Dl does not admit a 
decomposition of the form ® This will constitute a proof that the classical 

Fokker-Planck equation for the Lorentz plasma is not Hamiltonian. Suppose that does admit 
such a decomposition. Then the symmetric covariant tensor on phase space ul given by 

y) = iPx^o) ■ Dl ■ (iyi^o), 

■dhk){dhk-Y) 

k 


(5.134) 
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admits the reproducing kernel primitive k{zi,Z 2 ) = Yhk hk{z 2 ). Moreover, because aL can 
also be expressed as 


= i^{\v\)\Jijdx^dx^ (5.135) 

we also see that is positive semi-definite and has the constant rank 2 whenever |r;| 7 ^ 0. There¬ 
fore by Theorem 1, we must have Lxchl = 0 for any vector field X that takes values in the kernel 
of ai- On the other hand, one such X is given by X = Xh^ = v ■ d^, and 

(Lx„aMYiZ2) 

= L.v„. (oi(yi, I'j)) - ai .{ Y „ Lx„,l2) - aL ( Lx „ Yi , I2), (5.136) 

for arbitrary vector fields Yi, Y 2 on phase space. In particular, when Yi = Y 2 = w ■ dx + w ■ dy 
where w is a constant 3-component vector, 

{Lx,^aL){Y,, Y2 ) = 2wU- w, (5.137) 

which is never zero everywhere in phase space. This contradiction implies that ai does not admit 
a reproducing kernel primitive, which in turn implies that does not admit a decomposition of 
the form Dl = ® Xh ,.. Thus, the classical Fokker-Planck equation for the Lorentz plasma 

is not Hamiltonian. 

On energy-conserving Hamiltonian collision operators 

While Cl is not Hamiltonian and Chl does not conserve energy, perhaps there is some other 
collision operator C'(/) = drv{D-df) that approximates Cl, satisfies the Hamiltonian property, and 
conserves energy. The method introduced in BZQ would not be enough to find such an operator. 
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but it is still worthwhile to ask if such an operator exists, even in principle. In this section we 
will show that that the existence of such a C is unlikely in the following sense. We will say 
that the collision operator, div(ii) ■ df) in a Fokker-Planck equation conserves energy path-wise 
if the corresponding diffusion tensor satisfies D ■ dHo = 0, i.e. dHo is a null eigenvector of D. 
Note that Cl conserves energy path-wise. We will show that any Hamiltonian collision operator 
C that conserves energy path-wise is necessarily quite different from Cl, i-C- C must be a poor 
approximation to Cl- 

Let D = Yhk ^hk ® be the diffusion tensor associated with the Hamiltonian collision 
operator C that conserves energy path-wise. Associated with D is the symmetric covariant tensor 
a = Because D ■ dHo = 0, a ■ Xh^ = 0. In particular, 

0 = = Y,{hk,Ho}\ (5.138) 

k 

which implies that each of the Poisson commute with Ho- Therefore, 

Lxhq^ = 2 d{{hk, Hq}) dhk = 0. (5.139) 

k 

A reading of the proof of Theorem 3 now shows that there must be a symmetric covariant tensor 
d defined on the space of free streaming trajectories that pulls back to give a. A free streaming 
trajectory is a subset of TQ of the form 


l{xo,vo) = {{x, v) e TQ \ v = Vo and G M s.t. x = Vot + Xo}, (5.140) 

where {xo,Vo) is an arbitrary point in TQ with v 0. The space of free streaming trajectories, 
FS, is simply the union of all free streaming trajectories. Because the free streaming trajectories 
with a given velocity v can be identified with points in the plane perpendicular to v, FS has the 
structure of a rank-2 vector bundle over the 3-dimensional velocity space with the zero velocity 
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excluded, — {0}. To be precise, FS is diffeomorphic to the subbundle of ^ given 

by 


FS = {(n, x) G Mg X I a; • L = 0}. (5.141) 

There is also a natural projection map tt : TQ —)■ FS given by 7 r(a;, v) = (n, x±), where = 
X — X ■ vv. We must have a = tt* a, where d is some symmetric covariant tensor on FS. 

We will now analyze the form of d given that it must be the pullback of d along tt. Let 
ei(t^),e2(L) be orthogonal unit vectors that are everywhere perpendicular to v, i.e. ei{v) ■ v = 
e2(v) ■ V = 0. Using these unit vectors, we can introduce local coordinates on FS, (v,xi,X2), 
given by 


V = V (5.142) 

xi = X ■ ei(v) (5.143) 

X2 = x-e2(v), (5.144) 


where (v,x) G FS. If X = ■ d^+X^ • is a vector on the velocity phase space, its pushforward 

along TT is given by 


Tn{X) = • Cl + ^ ■ (xyci + |n|i?e2 • x) 

X"^ ■ 62 + • (a;||e2 - \v\Rei ■ x) 

d 


d 

dxi 

d 

dX2 


ov 


(5.145) 


where R = (V„ei) • 62 and x\\ = x ■ v. Now because 


d(Xi,X2) = d(T7r(Xi),T7r(X2)), 


(5.146) 
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and a is independent of x\\, we see that the eomponents of a must beeome large as x\\ beeomes 
large. Beeause Di does not depend on x, let alone x\\, by adjusting x\\, D can always be made 
mueh larger than D^. This rules out the possibility that D eould be a good approximation to Dl. 

5.2.6 Outlook 

We have suceeeded in showing that there is a Hamiltonian Fokker-Planck equation that governs the 
dynamics of the Lorentz plasma. This suggests that pitch angle scattering dynamics is Hamiltonian 
in a stochastic sense. On the other hand, we have also shown that our Hamiltonian Fo kk er-Planck 
equation does not have an exact energy conservation law. Instead the mean kinetic energy grows 
on a time scale proportional to r/rac, where r is the Fokker-Planck time step and t^c = (^p^- 

When inquiring as to why Chl produces energy, it is useful to remember the following fact 
about the derivation of Fokker-Planck equations. In these derivations, one tacitly (or explicitly, 
as we have done) divides the time axis into r-second intervals and studies the dynamics of an 
electron on a typical such interval. Because the positions of the ions are assumed to be statistically 
independent, this problem is reduced to studying the interaction of a single electron with a single 
Debye screened ion on a r-second time interval. Given the location of the ion, most electrons will 
both begin and end their r-second journeys well outside of the support of the ion’s potential. Each 
of these electrons will not suffer any change in its kinetic energy. However, there are some electrons 
that will either start or end within the support of the ion’s potential. Each of these electrons will 
suffer a change in their kinetic energy as a result of either climbing out of or falling into the ion’s 
potential well. Thus, an ensemble of electrons will not precisely conserve its kinetic energy over a 
r-second time interval. 

While this reasoning sheds some light on the energy-production problem, it is still not com¬ 
pletely satisfactory. Yes, an ensemble of electrons will only approximately conserve its kinetic 
energy over r-second time intervals. However, it is not physically true that the ensemble’s mean 
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kinetic energy can grow without bound, in contrast with the long-time behavior predicted by Chl 
(see Eq. (5.1 17)). This would correspond to electron heating, which is impossible because there is 
no energy impinging on the Lorentz plasma. 

The key to a complete understanding of the energetics of CHL is a careful scrutinization of the 
Markov assumption. This assumption artificially eliminates any memory an electron might have of 
its past after each r-second time interval. As a result, within the Markov model, the following non¬ 
physical dynamical process is possible. After moving for r seconds, an electron finds itself within 
the support of an ion’s potential. Whereas this electron should begin the next r-second interval by 
climbing out of this potential, instead it forgets the locations of all ions during the previous step, 
and, with high probability, fails to shed any of the kinetic energy it gained. This type of unphysical 
behavior allowed within the context of the Markov model is ultimately the source of the artificial 
heating predicted by Eq. (5.117). Moreover, it can be shown that this heating is not present in the 
classical Eokker-Planck equation for the Eorentz plasma because the probability that one of the 
unphysical processes just discussed occurs tends to zero as Tac/r tends to zero. 

We believe this explanation of why Chl causes slow artificial heating suggests a way to rec¬ 
oncile the stochastic Hamiltonian approach with energy conservation. The idea is to slightly relax 
the Markov assumption. In the rare event that electron ends a r-second time interval within the 
support of an ion’s potential, then the electrostatic field it sees in the next r-second interval should 
consist of the the old ion’s potential plus the potential of a new ion. Each time a new ion is drawn 
at random, it should never been chosen so that the support of its potential intersects the electron’s 
location. By allowing for this small memory effect, the unphysical electron trajectories that cause 
Chl to produce energy would be eliminated. Moreover, the evolution of a given electron would 
still be given by iterated symplectic mappings on phase space, and this ought to keep things within 
the realm of stochastic Hamiltonian mechanics. A challenging, yet enticing aspect of this possible 
route to overcoming the shortcomings of this Chapter would be developing the theory of stochastic 
Hamiltonian processes with memory. 
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